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with  the  objective  being  to  obtain  the  best  overall  fit  to  a 
given  experimental  relation  or  set  of  observed  responses.  As  a 
result,  the  accuracy  and  efficiency  of  this  model  calibration 
process  can  be  highly  dependent  on  both  the  subjectivity  of  the 
user  as  well  as  his  familiarity  and  expertise  with  the  particu¬ 
lar  constitutive  formulation.  ^In  order  to  ininimize-this  user-*’ 
dependence,  and  thereby  significantly  reduce'  the  complexity  of 
the  model  calibration  process,  a  computer-aided  automated 
procedure  has  been  developed  and  tested.  The  computer  code 
employs  a  Quasi-Newton  optimization  strategy  to  locate  that  set 
of  parameter  values  which  minimizes  the  discrepancy  between  the 
model  predictions  and  the  experimental  observations  included  in 
the  calibration  data  base.  Through  application  to  a  number  of 
real  soils,  the  automated  procedure  has  been  found  to  be  an 
efficient,  reliable  and  economical  means  of  accomplishing  model 
calibration.  Although  the  code  was  developed  specifically  for 
use  with  the  Bounding  Surface  plasticity  model,  it  can  readily 
be  adapted  to  other  constitutive  formulations.  Since  the  code 
greatly  reduces  the  dependence  of  calibration  success  on  user 
expertise,  it  significantly  increases  the  accessibility  and 
usefulness  of  sophisticated  material  models  to  the  general 
engineering  community 
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L  STRUCTURE  AND  CAPABILITIES 
1.1)  Introduction 

Much  research  has  been  directed  in  recent  years  towards  the  development 
of  sophisticated  constitutive  models  which  can  more  accurately  account  for  the 
diverse  stress>strain  phenomena  exhibited  by  soils  and  other  earth  materials. 
One  of  the  most  promising  of  these  new  models  is  a  plasticity  based  formulation 
known  as  the  Bounding  Surface  model.  The  Bounding  Surface  formulation  was 
originally  introduced  for  metals  by  Daf alias  and  Popov  (1974,  1975)  and  Dafalias 
(1975),  and  later  extended  to  cohesive  soils  by  Dafalias  (1979)  and  Dafalias  and 
Herrmann  (1980,  1982).  The  model  is  built  within  the  framework  of  traditional 
critical  state  soil  mechanics  and  employs  the  concept  of  a  bounding  surface  in 
stress  space.  It  has  been  shown  to  have  the  ability  to  accurately  capture  both 
the  drained  and  wdrained  behavior  of  clay  type  soils,  at  any  overconsolidation 
ratio,  under  either  monotonic  or  cyclic  loading. 

In  its  mast  general  form,  the  Boinding  Surface  model  requires  the 
determination  of  19  separate  constitutive  parameters,  including  2  initial  state 
properties,  5  traditional  material  constants,  whose  values  may  be  directly  obtained 
from  simple  well  known  laboratory  experiments,  and  12  model  constants,  which 
must  be  indirectly  established  through  a  trial  and  error  curve  fitting  process 
using  the  results  of  undrained  triaxial  testing.  A  general  summary  of  the  various 
properties  is  presented  by  Herrmann  et  al  (1980),  and  a  more  detailed  description 
of  both  the  qualitative  and  quantitative  influence  of  each  parameter  is  provided 
by  DeNatale  (1982). 

This  breakdown  of  model  constants  is  common  to  most,  if  not  all,  of  the 
soil  formulations  introduced  in  recent  years.  Determination  of  the  directly 
measureable  or  "fixed"  parameters  is  straightforward  and  readily  accomplished. 
Determination  of  the  remaining  "free"  parameters,  however,  can  make  the 
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calibration  procedure  prohibitively  difficult.  Rather  than  being  measured  directly 
from  a  particular  portion  of  a  specific  laboratory  test,  these  so-called  "free" 
parameters  must  be  established  by  trial  and  error,  with  the  objective  being  to 
obtain  the  best  overall  fit  to  a  given  experimental  relation  or  set  of  observed 
responses.  As  a  result,  the  accuracy  and  efficiency  of  the  calibration  process 
can  be  strongly  dependent  on  both  the  subjectivity  of  the  user  as  well  as  his 
expertise  with  the  particular  material  model. 

In  formulations  such  as  the  Bounding  Surface  model,  which  employ  a  small 
number  of  material  parameters  whose  roles  in  the  constitutive  formulation  are 
each  well  defined,  the  calibration  process  becomes  systematic  and 
straightforward.  However,  reliance  on  user  expertise  is  still  high,  since  all 
manual  curve  fitting  procedures,  by  their  very  nature,  require  both  judgement 
On  deciding  just  what  constitutes  the  "best"  overall  fit)  and  familiarity  On 
deciding  how  much  each  parameter's  value  must  be  changed  to  improve  a  given 
prediction). 

In  order  to  simplify  the  model  calibration  process,  a  computer  code  has 
recently  been  developed  by  DeNatale  (1982)  and  tested  on  a  variety  of  real 
soils.  The  code  employs  a  Quasi-Newton  optimization  strategy  to  locate  that 
set  of  parameter  values  which  minimizes  the  discrepancy  between  the  model 
predictions  and  the  experimental  observations  included  in  the  calibration  data 
base.  Because  this  new  computer  aided  procedure  greatly  reduces  the  dependence 
of  calibration  success  on  user  expertise,  it  significantly  increases  the  accessibility 
and  usefulness  of  sophisticated  material  models  to  the  general  engineering 
community. 

1.2)  The  Objective  Fwction 

Since  the  calibration  of  a  constitutive  model  involves  minimizing  the 
error,  or  residual,  between  the  observed  and  predicted  material  response,  the 


process  can  quite  naturally  be  viewed  as  an  optimization  problem.  Hence,  in 
‘order  to  develop  a  computer  directed  calibration  procedure,  it  is  necessary  to 
(i)  construct  an  objective  function  to  serve  as  a  scalar  measure  of  the  goodness 
of  a  particular  solution,  and  (ii)  select  a  search  strategy  to  enable  the  minimum 
of  this  function  to  be  located  in  an  efficient  and  reliable  manner.  In  both  the 
computer  code  and  the  discussion  to  follow,  it  is  assumed  that  the  model 
predictions  are  closely  enough  spaced  so  that,  by  joining  adjacent  points  with 
Unear  segments,  it  is  possible  to  obtain  a  good  approximation  to  a  smooth  curve 
(and  thus  the  concept  of  a  prediction  "curve"  is  vaUd). 

In  forming  an  objective  function,  the  total  error  between  the  experimental 
observations  and  model  predictions  could  be  expressed  by  either  (i)  summing  the 
discrete  residuals  at  each  of  the  experimental  points  included  in  the  calibration 
data  base,  or  (ii)  for  each  response  relation,  fitting  smooth  polynomials  through 
both  the  observed  and  predicted  data,  and  then  integrating  numerically  to  obtain 
the  area  (or  residual)  between  the  two  curves.  In  the  present  study  the  "discrete" 
approach  was  employed,  since  it  appeared  to  be  the  most  efficient  and  versatile 
of  the  two. 

In  most  data  fitting  and  regression  routines,  the  problem  is  generally 
posed  in  terms  of  one  or  more  independent  variables  x  and  a  dependent  variable 
y=y(x).  The  best  fit  to  the  data  is  then  obtained  by  minimizing  some  scalar 
measure  which  reflects  the  error  in  the  dependent  variable  y.  This  definition 
of  error  will  hereafter  be  referred  to  as  the  "vertical"  measure. 

In  geotechnical  engineering,  a  number  of  different  variables  are  typically 
used  to  define  the  soil  response.  In  a  conventional  imdrained  triaxial  compression 
test,  for  example,  the  response  might  be  characterized  in  terms  of  such  quantities 
as  the  mean  normal  effective  stress  p“,  the  deviator  stress  q,  the  pore  water 
stress  u,  and  the  axial  strain  ej*  The  distinction  between  independent  and 
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dependent  parameters  is  not  at  all  clear.  In  a  so-called  "strain-controlled**  test, 
it  could  perhaps  be  argued  that  represents  the  independent  variable;  therefore 
a  vertical  measure  may  be  appropriate  when  considering  such  relations  as  q  vs 
Ej  and  u  vs  Ej.  But  if  the  data  were  portrayed  as  a  stress  path  in  q-p"  space, 
both  axes  would  then  represent  dependent  parameters,  and  the  traditional 
classification  would  again  break  down. 

For  those  cases  in  which  all  aspects  of  the  soil  response  can  be  expressed 
in  terms  of  a  single  quantity  (by  using,  for  example,  p'  vs  Ej,  q  vs  Ej  and  u 
vs  Ej  rather  than  q  vs  p',  q  vs  Ej  and  u  vs  E|),  a  vertical  measure  may  be 
reasonable.  However,  if  the  distinction  between  independent  and  dependent 
variables  cannot  be  made  clear,  there  is  no  more  reason  to  use  a  vertical 
measure  (y=y(x))  than  there  is  to  use  a  horizontal  measure  (x=x(y)).  In  these 
cases  it  is  probably  more  appropriate  to  use  a  measure  such  as  the  shortest 
distance  between  the  experimental  observation  and  prediction  curve.  While  there 
is  no  reason  to  suspect  that  this  alternative  is  theoretically  more  sound,  such 
a  "Euclidean"  measure  is  probably  closer  to  what  one  intuitively  uses  when 
estimating,  by  sight,  the  error  between  two  curves. 

In  addition  to  defining  a  direction  along  which  the  error  will  be  measured, 
it  is  also  necessary  to  consider  whether  the  sum  of  the  signed,  absolute  or 
squared  residuals  should  be  minimized.  In  regression  analysis,  a  squared  measure 
has  traditionally  been  used.  The  squared  measure  not  only  yields  a  unique 
solution,  but  also  has  additional  statistical  relevance  (as  described,  for  example, 
by  Alder  and  Roessler  (1968)).  The  absolute  measure  likewise  yields  a  unique 
solution,  but  is  generally  avoided  due  to  the  mathematical  difficulties  customarily 
associated  with  the  occurance  of  absolute  values.  The  signed  measure  is  only 
rarely  used,  because  it  can  lead  to  non-unique  solutions. 

Since  there  is  no  theoretical  reason  to  select  one  or  the  other,  the 
calibration  code  permits  ei*  >cr  absoiut*  or  squared  residuals  to  be  used  to  form 
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the  objective  function.  When  these  two  options  are  combined  with  the  choice 
of  either  vertical  or  Euclidean  measure  (as  described  above),  the  user  then  has 
four  possible  ways  to  define  the  error  at  a  point.  Recent  research  by  DeNatale 
(1982)  has  shown  that  the  location  of  the  global  minimum  remains  essentially 
the  same,  regardless  of  which  options  are  used.  However,  preliminary  applications 
to  a  variety  of  artificial  and  real  soils  indicate  that  the  absolute-Euclidean 
measure  results  in  a  more  well-behaved  objective  function  that  can  be  most 
easily  minimized. 

Mi: 

In  order  to  combine  the  residuals  at  various  points  within  a  given  relation 
(for  example,  the  q  vs  Cj  relation  from  a  test  at  OCR=l),  or  from  relations  of 
tWfc'  same  kind  obtained  from  different  tests  (for  example,  the  u  vs  Cj  relation 
from  tests  at  OCR=l  and  2),  it  is  necessary  to  first  define  what  is  meant  by 
"equal"  error.  A  given  solution  could  be  defined  to  be  equally  good  at  two 
points  A  and  B  if  either  (i)  there  was  the  same  absolute  r*ror  at  both  A  and 
B,  or  (ii)  there  was  the  same  relative  error  at  A  as  at  B.  The  decision  as  to 
which  definition  of  equality  should  be  used  is  entirely  up  to  the  user.  The 
computer  code  permits  the  use  of  either  extreme,  or  any  point  in  between. 
Consequences  of  the  various  choices  are  discussed  further  by  DeNatale  (1982). 
It  should  be  noted  that  most  physical  analogies  lie  somewhere  in  the  middle  —a 
given  dial  gauge  or  pore  pressure  transducer  may  be  accurate  to  within  a%,  but 
may  fail  to  register  any  meaningful  readings  below  a  magnitude  of  8. 

The  code  permits  any  number  of  tests,  relations  and/or  individual 
observations  to  be  included  in  the  calibration  data  base.  Because  the  various 
response  relations  will  generally  be  of  different  dimensions  (such  as  stress  vs 
stress,  stress  vs  strain,  strain  vs  strain,  etc.),  all  data  is  nondimensionalized  so 
that  errors  from  different  relations  can  be  properly  combined.  Different  weights 
may  be  assigned  to  specific  components  of  the  data  base  if  it  is  felt  that  certain 


tests,  relations  or  observations  are  more  reliable  or  representative  than  others, 
or  if  it  is  necessary  to  have  the  final  model  predictions  fit  some  data  more 
closely  than  others.  T.\a  consequences  and  proper  role  of  weighting  factors  is 
again  discussed  by  DeNatale  (1982). 

1.3)  The  Optimization  Strategy 

An  extremely  large  number  of  optimization  strategies  have  been  suggested 
over  the  last  30  years,  with  the  performance  of  a  given  approach  being  strongly 
dependent  on  the  particular  type  of  problem  to  which  it  is  applied.  Hence,  it 
is  generally  agreed  that  there  is  no  single  best  algorithm,  but,  rather,  only 
strategies  which  perform  best  when  applied  to  certain  classes  of  problems. 

In  selecting  the  most  suitable  approach,  a  key  factor  is  whether  or  not 
first  and  second  derivative  information  can  be  readily  obtained.  With  the 
Bounding  Surface  model,  as  with  most  sophisticated  material  models,  the 
governing  equations  are  so  complex  that  it  is  essentially  impossible  to  directly 
relate  the  relevant  response  parameters  (such  as  p ",  q,  u,  e},  etc.)  to  the 
constitutive  parameters  employed  by  the  formulation.  Thus,  the  objective 
function  must  be  formed  by  sunming  a  series  of  discrete  weighted  residuals, 
and  therefore  first  and  second  derivative  information  is  not  explicitly  available. 

Under  these  conditions,  the  current  consensus  among  those  most  active 
in  optimization  research  (including,  for  example,  Fletcher  (1980)  and  Gill  et  al 
(1981))  is  that  a  Quasi-Newton  strategy  with  finite  difference  approximations  to 
derivatives  will,  if  properly  implemented,  generally  exhibit  the  most  efficient 
and  reliable  performance.  Hence,  a  Quasi-Newton  strategy  was  incorporated 
into  the  calibration  code  to  direct  the  search. 

The  Quasi -Newton  algorithms  avoid  most  of  the  difficulties  associated 
with  Newton’s  method  (as  outlined,  for  example,  by  Murray  (1972)),  while  retaining 
the  roubstness  of  the  second  derivative  methods.  In  contrast  to  Newton's  method, 
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where  the  Hessian  G  is  evaluated  directly  at  e^ch  iteration,  the  Quasi-Newton 
routines  build  up  curvature  information  gradually  as  the  search  proceeds,  relying 
only  on  the  observed  behavior  of  the  objective  function  f  and  its  gradient  g. 

In  the  Quasi-Newton  approach,  the  inverse  Hessian  G~*  is  approximated 
by  a  symmetric,  positive-definite  matrix  H  which  is  updated  at  each  iteration. 


The  kth  iteration  of  the  search  has  the  form  (Fletcher  (1980)): 


i)  set  sK 


-Hkgk 


K 

ii)  perform  a  line  search  along  s  ,  ending  at 

k+1  k  k  k 
xK+1  =  xK  +  as 

iii)  update  Hk,  yielding  Hk+* 

k  k 

where  x  denotes  the  kth  solution  vector,  s  represents  the  kth  search  direction, 
and  a  is  a  positive  scalar. 

Since  the  basic  premise  in  the  Quasi-Newton  strategy  is  that  curvature 
information  can  be  estimated  without  explicitly  forming  the  Hessian  G,  most 
research  with  these  methods  has  focused  on  the  development  of  improved  formulae 
for  updating  H  in  step  (iii)  above.  At  the  present,  the  BFGS 
(Broyden-Fletcher-Goldfarb-Shanno)  update  is  generally  regarded  as  the  most 
effective.  The  update  has  the  form: 

^  .  H  .  <1  .£*)<«£)  -  («TTH  *  mL ) 

Brcs  6  y 


where: 


k  k+1  k 

y  =  g  -  g 

-k  k+l  k  k  k 

6  =  x  -  x  =  a  s 


and  where  the  superscript  k  has  been  dropped  for  convenience.  Further  details 
of  the  BFGS  update  and  other  aspects  of  the  Quasi-Newton  methods  are  provided, 
for  example,  by  Broyden  (1972),  Dennis  and  Mori  (1977)  and  Fletcher  (1980). 


To  accomplish  the  line  search  in  step  (ii)  above,  a  strategy  outlined  by 
Fletcher  (1980)  was  employed,  together  with  certain  ad  hoc  modifications 
described  by  DeNatale  (1982).  This  particular  strategy  permits  the  line  search 
to  be  carried  out  to  any  degree  of  accuracy,  which  is  ideal,  since  a  primary 
advantage  of  the  BFGS  update  is  that  it  generally  performs  most  efficiently 
with  low  accuracy  line  searches. 

Because  no  analytic  expressions  are  available  for  first  derivatives  g, 
gradients  must  be  evaluated  by  finite  differences.  The  strategy  used  is  patterned 
after  that  described  by  Stewart  (1967).  Directed  by  the  local  geometry  of  the 
objective  function,  the  algorithm  switches  between  forward  and  central 
differencing  formulae,  and  continually  adjusts  the  differencing  interval,  in  an 
attempt  to  balance  round-off  and  truncation  errors. 

Finally,  the  code  permits  constraints  to  be  imposed  on  the  various 
constitutive  parameters  in  the  form  of  simple  bounds: 

4i  <  xi  <  ui 

where  JL  and  Uj  represent  the  minimum  and  maximum  values  that  parameter  x. 
can  assume.  A  restriction  of  this  kind  would  be  appropriate  if  either  0)  there 
were  certain  theoretical  restrictions  placed  on  the  value  which  a  given  parameter 
could  assume,  or  (ii)  there  were  certain  ranges  of  parameter  values  beyond  which 
the  numerical  implementation  of  the  material  model  became  unstable,  or  (iii) 
certain  material  properties  were  observed  experimentally  to  vary  over  some 
finite  range,  and  there  was  no  overwhelming  reason  to  fix  them  prior  to  calibration 
at  any  particular  values. 

There  are  a  number  of  ways  to  directly  (through  the  use  of  a  constrained 
optimization  code)  or  indirectly  (through  the  use  of  barrier  and  penalty  functions, 
Lagrange  multipliers  and/or  variable  transformations)  account  for  simple  bounds. 


After  a  review  of  available  literature,  it  was  concluded  that  a  trigonometric 
variable  transformation  would  provide  the  most  straightforward  solution.  Hence, 
rather  than  minimizing  f(x)  subject  to  t.<x.<u.,  the  program  makes  the 
transformation: 

*i  *  *i  ♦  sin2  y. 

or: 

yA  r  sin"1  [(xrli)/(urg.)l1/2 

and  minimizes  f  [x(y)]=f(y),  where  y.  can  now  assume  any  value  in  the  range 
-«*><y.  <+<*>.  Although  variable  transformations  are  potentially  subject  to  a  number 
of  difficulties  (as  described,  for  example,  by  Gill  et  al  (1981)),  the  above 
transformation  was  found  to  work  extremely  well  in  this  application,  with  no 
detectable  problems. 

!.%)  Generation  of  Model  Predictions 

Tne  search  for  the  optimal  set  of  parameter  values  is  directed  by  the 
Quasi-Newton  strategy  previously  described.  However,  in  order  to  evaluate  the 
objective  function  at  some  location  x,  it  is  first  necessary  to  generate  a 
corresponding  set  of  model  predictions.  To  accomplish  this,  the  calibration  code 
relies  on  two  subroutines  —  EVAL  and  CLAY  —  developed  by  Herrmann  et  al 
(1980,  1981)  during  previous  research  with  the  Bounding  Surface  model.  Subroutine 
EVAL  performs,  essentially,  single-element,  incremental-iterative  finite  element 
analyses  of  bodies  under  a  homogeneous  state  of  stress  and/or  strain.  Subroutine 
CLAY  consists  of  a  numerical  implementation  of  the  governing  constitutive 
equations,  and  thus,  when  called,  provides  the  appropriate  material  response  to 
the  given  stress  and/or  strain  increment.  In  order  to  adapt  the  present  calibration 
code  to  other  constitutive  models,  one  need  only  replace  CLAY  with  the 
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appropriate  materials  subroutine.  Additional  details  of  subroutines  EVAL  and 
CLAY  are  presented  by  Herrmann  et  al  (1980,  1981,  1982). 

1.5)  The  Calibration  Data  Base 

The  ultimate  goal  of  the  calibration  process  is  to  identify  that  set  of 
parameter  values  which  enables  the  theoretical  model  to  most  closely  simulate 
the  observed  material  response.  This  goal  is  ordinarily  accomplished  by  fitting 
the  model  to  a  representative  set  of  experimentally  observed  stress-strain 
relations  or  "calibration  data  base."  Ideally,  this  calibration  data  base  should 
be  complete  and  diverse  enough  that  all  important  aspects  of  the  material's 
response  are  included,  and  all  necessary  constitutive  parameters  may  be  uniquely 
established. 

In  its  most  general  form,  the  Bounding  Surface  formulation  becomes  a 
fully  three-dimensional  stress-strain  modeL  With  a  single  set  of  parameter 
values,  the  model  may  be  applied  to  specimens  at  all  overconsolidation  ratios 
(OCR's),  subjected  to  either  monotonic  or  cyclic  compression  and/or  extension 
loading,  mder  either  drained  or  tndrained  conditions.  Hence,  to  establish  the 
optimal  values  of  the  necessary  constitutive  parameters,  the  calibration  data 
base  should  ideally  contain  observations  (in  the  form  of  q  vs  p',  q  vs  tj  and 
u  vs  G|  relations)  from  the  following  seven  standard  laboratory  tests: 

i)  isotropic  (or  Kq)  consolidation  or  drained  compression 

test,  with  both  loading  and  unloading;  and, 
ii-vii)  undrained  triaxial  compression  and  extension  tests  on 
specimens  in  the  normally  (OCR  =  1),  lightly  (1  <  OCR  <  2.5) 
and  heavily  (OCR  >  4)  over  consolidated  regions. 

The  results  of  the  consolidation  test  are  required  to  establish  the  slopes  of  the 
isotropic  consolidation  and  swell /recompression  curves  in  e-In  p"  space,  X  and  k. 
These  two  parameters  belong  to  the  class  of  traditional  material  properties,  and 


would  normally  be  assigned  values  immediately,  prior  to  using  the  automated 
calibration  procedure.  The  results  of  the  six  undrained  triaxial  experiments  are 
required  to  determine  the  12  model  constants  cited  in  section  1.1,  and  would 
thus  provide  the  data  needed  to  direct  the  automated  calibration  procedure. 
Naturally,  if  a  less  general  form  of  the  Bounding  Surface  model  is  acceptable, 
the  number  of  constitutive  parameters  involved,  and  the  number  of  laboratory 
experiments  required,  can  be  (tastically  reduced.  For  example,  if  the  model  is 
only  to  be  applied  to  normally  consolidated  soils  loaded  in  triaxial  compression, 
the  number  of  required  constitutive  parameters  (taps  from  19  to  7,  and  only 
the  isotropic  consolidation  results  and  a  single  triaxial  test  are  needed  for  model 
calibration. 

Although  the  above  data  base  is  strongly  recommended,  the  Bounding 
Surface  model  could  also  be  calibrated  using  other  types  of  data.  For  example, 
drained  rather  tht.y,  ndrained  tests  could  be  employed.  However,  imdrained 
tests  are  preferable,  since  good  initial  estimates  for  many  of  the  model 
parameters  can  be  made  by  examining  the  experimentally  observed  undrained 
stress  paths. 

There  is  also  some  evidence  that  the  calibration  data  base  need  not 
necessarily  include  data  from  all  three  regions  of  overconsolidation  ratio  (see 
DeNatale  (1982)).  That  is,  it  may  be  sufficient  to  include  only  tests  from  the 
normal  and  heavy  ranges,  or,  perhaps,  even  from  the  heavy  range  alone.  The 
data  which  supports  this  possibility  is  not,  however,  conclusive,  and  therefore 
testing  at  all  three  overconsolidation  regions  is  still  advised. 

In  addition,  the  experimental  observations  need  not  necessarily  include  all 
three  response  relations  q  vs  p",  q  vs  Ej  and  u  vs  Cj.  Note  that  of  the  four 
undrained  response  parameters  p',  q,  u  and  £j,  only  three  are  independent.  In 
practice,  p"  is  never  actually  measured,  but,  rather,  is  computed  from  the 
relation: 
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where  represents  the  applied  lateral  confining  pressure.  Thus,  any  two  of 
the  three  relations  cited  above  will  completely  define  the  soil  response.  The 
use  of  q  vs  p'  or  q  vs  data  alone  is  insufficient,  since  each  of  these  relations 
is  insensitive  to  certain  of  the  constitutive  parameters.  There  is  some  evidence 
that  the  use  of  u  vs  data  alone  may,  however,  be  adequate  (see  DeNatale 
(1982)).  Nevertheless,  the  use  of  all  three  response  relations  appears  to  increase 
the  rapidity  with  which  the  optimization  algorithm  converges  to  the  minimum. 
Presumably,  the  inclusion  of  redundant  data  reinforces  the  correct  search 
direction.  Hence,  since  the  cost  of  a  computer  directed  calibration  run  is  only 
marginally  affected  by  the  number  of  response  relations  included  in  the  calibration 
data  base,  it  is  therefore  still  recommended  that  all  three  of  the  relations  cited 
above  be  used. 

Finally,  it  may  be  passible  to  use  testing  devices  other  than  the  triaxial 
apparatus  to  acquire  the  necessary  experimental  observations.  Although  the 
conventional  triaxial  apparatus  is  the  most  common  and  versatile  laboratory 
device,  the  simple  shear  apparatus  could,  for  example,  also  be  used.  In  general, 
the  material's  observed  stress-strain  characteristics  will,  to  some  extent,  be 
dependent  on  the  testing  device  employed.  Thus,  in  practical  problems,  the 
laboratory  device  used  to  acquire  the  calibration  data  base  should  simulate,  as 
closely  as  possible,  the  loading  conditions  for  which  Bounding  Surface  predictions 
will  eventually  be  generated. 

1.6)  Practical  Considerations 

As  previously  mentioned  in  section  1.2,  the  automated  calibration  code 
seeks  to  locate  the  objective  function's  global  minimum.  Unfortunately,  there 
is  no  guarantee  that  the  algorithm  will  always  succeed.  The  Quasi-Newton 


strategy  employed  by  the  model  calibration  code,  like  most,  if  not  all,  practical 
optimization  algorithms,  is  designed  only  to  locate  local  minima  in  the  vicinity 
of  the  initial  estimates.  Hence,  the  probability  that  the  true  global  minimum 
will  be  found  is  directly  related  to  the  degree  of  unimodality  exhibited  by  the 
objective  function  and  the  accuracy  of  the  initial  starting  guess. 

Preliminary  research  by  DeNatale  (19S2)  has  shown  that  the  use  of  the 
absolute-Euclidean  measure  of  error  leads  to  a  more  ixumodal,  and  thus  desireable, 
objective  function.  A  procedure  for  acquiring  improved  starting  estimates  has 
also  been  developed  by  DeNatale  (1982).  Through  actual  testing  with  a  number 
of  different  soils,  this  strategy  has  been  found  to  produce  starting  estimates 
which  enable  the  automated  calibration  code  to  consistently  locate  the  global 
minimum.  In  practice,  however,  the  only  way  to  ensure  that  the  global  minimum 
has  been  found  is  to  conduct  the  search  from  a  variety  of  different  starting 
points.  The  solution  which  yields  the  lowest  value  of  the  objective  function 
may  then  be  regarded  as  the  true  global  minimum. 

A  second  practical  consideration  concerns  the  quality  of  the  calibration 
data  base.  The  user  should  ensure  that  the  experimental  observations  included 
in  the  calibration  data  base  are  diverse  enough  to  permit  the  optimal  values  of 
the  required  unknown  model  parameters  to  be  uniquely  defined.  For  example, 
if  the  code  is  used  to  identify  those  model  parameters  associated  with  the 
heavily  over  consolidated  material  response,  the  calibration  data  base  must  include 
observations  made  on  heavily  overconsoli dated  soil  specimens.  If  the  necessary 
experimental  data  is  not  included,  the  program  will  continue  to  execute,  but 
the  final  computed  "optimal"  values  of  the  undefined  parameters  will  be  very 
close  to  the  initial  estimates. 

The  major  consequence  of  an  inadequate  or  incomplete  calibration  data 
base  is  therefore  related  to  the  cost  of  the  analysis.  Certain  computational 


costs  increase  in  proportion  to  n^  (where  n  represents  the  number  of  parameters 
whose  optimal  values  are  being  sought),  and  a  single  gradient  evaluation  requires 
either  n  or  2n  additional  objective  function  evaluations,  depending  on  whether 
forward  or  central  differencing  formulae  are  used.  Thus,  to  minimize  the  cost 
of  the  analysis,  the  user  should  seek  to  identify  only  those  parameters  whose 
optimal  values  can  be  defined,  given  the  particular  experimental  data  base.  A 
comprehensive  discussion  of  the  influence  of  each  of  the  19  model  parameters 
is  provided  by  DeNatale  (19S2),  which  may  be  referred  to  if  any  uncertainty 
exists. 

1.7)  Verification 

In  order  to  verify  the  viability  of  the  new  computer  aided  calibration 
procedure,  the  method  was  applied  to  a  number  of  representative  data  bases, 
both  artificial  and  real.  The  outcome  of  these  studies  is  discussed  in  detail  by 
DeNatale  (1982).  Among  the  real  data  bases  to  which  the  automated  process 
was  applied  was  the  experimental  research  on  Kaolin  reported  by  Wroth  and 
Loudon  (1967).  The  Bounding  Surface  model  was  calibrated  on  the  basis  of 
uidrained  triaxial  compression  tests  on  samples  at  overconsolidation  ratios  of 
OCR  =  1.0,  1.5  and  4.5.  With  the  necessary  constitutive  parameters  having  thus 
been  fixed  at  their  optimal  values,  predictions  were  then  generated  for  (i) 
undrained  triaxial  compression  tests  at  OCR  =  1.2,  1.8,  2.5,  3.0  and  6.5,  (ii)  an 
mdrained  cyclic  triaxial  compression  test  at  OCR  =  1.0,  and  (iii)  a  drained 
triaxial  compression  test  at  OCR  =  1.3.  The  experimental  observations  and 
model  predictions  are  compared  in  Figures  1.1  through  1.3.  Further  details  of 
the  calibration  procedure  and  predictions  for  this  particular  data  base  are  provided 
by  both  DeNatale  (1982)  and  Herrmann  et  al  (1982), 
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:  Experimental  Observations  and  Model 
Predictions  For  the  Undrained  Triaxial 
Compression  Tests 
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Figure  1.3:  Experimental  Observations  And  Model 
Predictions  For  The  Drained  Triaxial 
Compression  Test  (OCR  =  1.3) 


3 


1.8)  Cost 

The  automated  calibration  code  has  been  written  in  FORTRAN  and 
implemented  on  both  an  LSI- 11/23  minicomputer,  as  well  as  a  VAX-11/780 
super -minicomputer.  The  cost  of  a  given  analysis  is  controlled  primarily  by  the 
number  of  distinct  experimental  tests  included  in  the  calibration  data  base  and 
the  number  of  constitutive  parameters  whose  optimal  values  are  being  sought. 
A  typical  computer  calibration,  such  as  that  reported  in  section  1.7  for  the  data 
of  Wroth  and  Loudon  (1967),  requires  from  200-400  objective  function  evaluations, 
or  about  30-60  minutes  of  VAX  CPU  time,  at  a  cost  of  approximately  $25.00- 
$50.00.  When  compared  to  the  expense  associated  with  the  acquisition  of  the 
experimental  calibration  data  base,  or  the  expense  of  subsequent  finite  element 
analyses  involving  the  calibrated  model,  the  cost  of  a  typical  computer  aided 
calibration  becomes  relatively  low.  Note  also  that  the  material  model  need 
only  be  calibrated  once  for  any  particular  soil,  regardless  of  the  variety  and 
number  of  finite  element  analyses  that  may  sii>sequently  be  performed. 

1.9)  Conclusion 

In  order  to  facilitate  the  calibration  of  sophisticated  constitutive  models, 
an  automated  FORTRAN  code  has  been  developed  and  tested.  The  code  employs 
a  Quasi-Newton  optimization  strategy  to  locate  that  set  of  parameter  values 
which  minimizes  the  weighted  residual  between  the  model  predictions  and  the 
experimental  calibration  data  base.  Through  application  to  a  number  of  real 
soils,  this  new  procedure  has  been  found  to  be  an  efficient,  reliable  and  economical 
means  of  accomplishing  model  calibration. 

Although  the  code  was  developed  specifically  for  use  with  the  Bounding 
Surface  plasticity  model,  it  can  be  readily  adapted  to  other  constitutive 
formulations.  Because  the  code  greatly  reduces  the  dependence  of  calibration 
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success  on  user  expertise,  it  significantly  increases  the  accessibility  and  usefulness 
of  sophisticated  material  models  to  the  general  engineering  community. 

In  addition  to  being  a  useful  practical  tool,  the  code  can  also  be  of  great 
value  to  the  theoretician  involved  in  model  development.  Since  use  of  the  code 
enables  various  sets  of  model  predictions  to  be  quantitatively  compared,  it 
enables  one  to  investigate  such  features  as  (i)  the  model's  sensitivity  to  each 
of  the  constitutive  parameters,  Gi)  the  significance  and  uniqueness  of  each 
parameter's  effect,  (iii)  the  influence  of  the  calibration  data  base  on  the  computed 
optimal  solution  and  Gv)  the  uniqueness  of  a  given  solution  or  set  of  model 
predictions.  These  uses  of  the  code,  as  well  as  several  others,  are  discussed 
and  illustrated  in  greater  detail  by  DeNatale  (1982). 

The  computer  aided  procedure  is  not,  however,  a  cure-all,  and  use  of  the 
automated  scheme  does  not  guarantee  that  all  problems  associated  with  manual 
model  calibration  will  disappear.  Many  of  the  difficulties  encountered  during 
the  manual  calibration  of  a  particular  material  model  have  little  or  nothing  to 
do  with  user  expertise,  but,  rather,  arise  due  to  the  nature  of  the  formulation 
itself.  The  presence,  in  a  particular  model,  of  an  excessive  number  of  constitutive 
parameters,  or  ill-defined  parameters  with  non-uiique  effects,  will  always  make 
a  unique  solution  difficult,  or  perhaps  impossible,  to  find,  regardless  of  whether 
a  computer  approach  is  used  or  not.  Thus,  it  is  best  to  regard  the  automated 
procedure  as  a  useful  tool  for  well-constructed  models,  rather  than  as  a  fool-proof 
cure  for  ill-constructed  ones. 
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INPUT  INSTRUCTIONS  FOR  MODCAL 


Heading  Information  Cards 
1st  Card  (20A»)s 


Columns 

1-80  TITL  :  any  information  that  is  to  be  printed  as 

a  title  for  the  analysis 

Control  Codes  and  Analytical  Options  Card: 

1st  Card  (1615): 


Columns 


1  -  5 

NTST 

• 

• 

number  of  distinct  experimental  tests  for 
which  model  predictions  are  desired,  or 
on  which  model  calibration  is  to  be  based 

6-10 

JOPT 

= 

0  -*•  model  calibration  not  required 

1  model  calibration  required 

11  -  15 

3RUN 

= 

0  ■*>  model  predictions  not  generated 

1  model  predictions  generated 

16  -  20 

JPLT 

0  observed  and  predicted  responses 

not  plotted 

1  ■*>  observed  and  predicted  responses 
plotted 

21  -  25 

KRSL 

= 

0  -*•  Euclidean  measure  of  error  used 

1  -*■  vertical  measure  of  error  used 

26  -  30 

KNRM 

s 

0  -*■  absolute  residuals  used 

1  squared  residuals  used 

31  -  35 

NOPT 

• 

• 

number  of  parameters  to  be  established 
by  optimization 

36-80  KOPT  (ID,  II  =  1,  NOPT:  property  numbers  of  the 

NOPT  parameters  to  be  established  by 
optimization  (see  section  Ill  and  Table 
2.1).  Use  15  formats  and  continue  on  a 
2nd  card  if  necessary  (that  is,  if 
NOPT  >  9) 


HL 


Material  and  Model  Properties  Cards:^ 
1st  Card  (8EI0.3): 


Columns 


1  -  10 

X 

11  -  20 

K 

21-30 

Mc 

31  -  40 

Me/M 

41-50 

G 

51  -  60 

r 

61  -  70  pt 


71  -  80  p, 

O 


slope  of  the  Isotropic  consolidation  line 
in  e  -  in  p  space 

slope  of  the  swell/recompression  line  in 
e  -  In  p  space 

slope  of  the  critical  state  line  in  triaxial 
compression 

ratio  of  extension  to  compression  values 

elastic  shear  modulus  (or,  alternatively, 
Poisson's  ratio  v) 

combined  bulk  modulus  of  the  pore  water 
and  soil  skeleton.  For  drained  conditions, 
T  =  0.  For  undrained  conditions,  .a  value 
in  the  range  T  =  (1(T  to  10b)*p  is 
recommended  a 

transitional  stress  at  which  the  drained 
compression  response  changes  from  linear 
in  e  -  p  space  to  linear  in  e  -  In  p  space. 
If  no  data  in  this  region  is  available,  a 
value  in  the  range  p,  =  (0.3  to  1.0)*p 
is  recommended  a 

atmospheric  pressure 


2nd  Card  (7E10.3)i 
Columns 
1-10  R 

c 

11  -  20  A„ 

c 


surface  shape  parameter  associated  with 
ellipse  1  in  compression 

surface  shape  parameter  associated  with 
the  hyperbola  in  compression 


^  If  model  calibration  is  required,  it  is  not  necessary  to  enter  estimates 
for  those  parameters  whose  optimal  values  are  being  sought. 
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21  -  30 

T 

• 

• 

surface  shape  parameter  associated  with 
ellipse  2  in  compression  and  extension. 
A  value  in  the  range  0.05  <  T  <  0.15 
is  suggested 

31  -  00 

V*c 

• 

• 

ratio  of  extension  to  compression  values 

01-50 

VAc 

m 

• 

ratio  of  extension  to  compression  values 

51  -  60 

c 

• 

• 

projection  point  parameter.  A  value  of 
c  =  0.0  places  the  projection  point  at 
the  origin 

61  -  70 

s 

• 

• 

elastic  nucleus  parameter.  A  value  of 
s  =  1.0  causes  the  elastic  zone  to 
degenerate  to  a  point 

3rd  Card  (0E1Q.3): 

Columns 

1  -  10 

he 

• 

• 

primary  hardening  parameter  in  compres¬ 
sion 

11  -  20 

mc 

• 

• 

secondary  hardening  parameter  in  com¬ 
pression.  A  value  in  the  range 

0.2<m<1.0  is  suggested 

21  -  30 

"A 

• 

• 

ratio  of  extension  to  compression  values 

31  -  00 

"Ac 

• 

• 

ratio  of  extension  to  compression  values 

»th  Cards  (15.  «,  4E103): 

One  card  is  required  for  each  of  the  NTST  distinct  experimental 
tests  for  which  model  predictions  are  desired: 

Columns 


1  -  5 

rrsT 

• 

• 

test  number 

11  -  20 

eo 

• 

• 

initial  void  ratio 

21  -  30 

Po 

• 

• 

initial  maximum  past  effective  precon¬ 
solidation  pressure 

o 

1 

°co 

• 

• 

initial  confining  pressure 

41-50 

OCR 

• 

• 

initial  over  consolidation  ratio,  as  defined 
by  OCR  =  po/oco 

Convergence  Criteria  and  Iteration  Information  Card: 


1st  Card  (315,  5X,  4E10.3): 
Columns 
1  -  5  KIND 

6-10  ITMX 

11-15  LARG 

21  -  30  CNFR 


31  -  40  ERMX 


41  -  50  PLIM 


51  -  60  TLIM 


0  ■*  reformulated  nearly  incompres¬ 
sible  analysis 

1  -*•  non-reform  ulated  nearly  incom¬ 

pressible  analysis 

maximum  number  of  iterations  permitted 
for  a  given  solution  increment.  Values 
in  the  range  5  <  ITMX  <  10  are 

suggested 

0  •*  engineering  stresses  and  strains 
assumed 

1  •*  true  stresses  and  natural  (log¬ 
arithmic)  strains  assumed 

establishes  CNFR  and  1/CNFR  as  the 
lower  and  upper  limits  for  the  calculated 
values  of  the  Aitken's  acceleration 
factors.  Values  in  the  range  0.0<  CNFR 

<  1.0  are  permitted,  and  if  CNFR  =  1.0, 
acceleration  factors  will  not  be  used 

maximum  allowable  relative  difference 
between  the  norms  of  the  incremental 
stress  and  strain  vectors  from  two  con¬ 
secutive  iterations.  Values  in  the  range 
0.001  <  ERMX  <  0.01  are  recommended. 
If  convergence  does  not  occur  within 
ITMX  iterations,  a  message  is  printed 
and  program  execution  will  cease 

percentage  of  the  maximum  absolute 
experimentally  observed  y-value  for  a 
given  plot,  below  which  absolute  rather 
than  relative  errors  are  used.  Values  in 
the  range  0.0  <  PLIM  <  1.0  are 

permitted,  with  PLIM  =  0.0  resulting  in 
relative  errors  being  used  at  all  points 
and  PLIM  =  1.0  resulting  in  absolute 
errors  being  used  at  all  points 

percentage  of  the  maximum  absolute 
experimentally  observed  y-value  for  a 
given  type  of  relation,  below  which 
absolute  rather  than  relative  errors  are 
used.  Values  in  the  range  0.0  <  TLIM 

<  1.0  are  permitted,  with  TLIM  =  0.0 
resulting  in  relative  errors  being  used  for 
all  plots  and  TLIM  =  1.0  resulting  in 
absolute  errors  being  used  for  all  plots 
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Specified  Loading  History  Cards: 


The  following  sequence  of  cards  must  be  supplied  for  each  of  the 
NTST  distinct  experimental  tests  for  which  model  predictions  are 
required: 

1st  Card  (215): 

Columns 

1-5  IT5T  :  test  number 

6-10  NSEC  :  number  of  distinct  history  segments  into 

which  the  test  is  subdivided 

2nd  Cards  (6  (11.  E9.2).  5X.  15.  E10.3): 

One  card  is  required  for  each  of  the  NSEG  distinct  history  segments 
into  which  the  test  is  subdivided: 


Columns 

1 

LTYP1 

= 

0  ♦  a  is  specified 

1  ♦  e*  is  specified 

2-10 

VALUj 

• 

• 

value  of  a  or  e  at  the  end  of  the  given 
history  segment1 

11 

ltyp2 

= 

0  ♦  o  is  specified 

1  -►  cy  is  specified 

12  -  20 

valu2 

• 

• 

value  of  o  or  c  at  the  end  of  the  given 
history  segment 

21 

LTYPj 

= 

0  ♦  o  is  specified 

1  +  e*  is  specified 

22  -  30 

VALUj 

« 

• 

value  of  o  or  c  at  the  end  of  the  given 
history  segment 

31 

LTYP# 

= 

0  -*■  t  is  specified 

1  ♦  is  specified 

xy 

© 

•a- 

i 

<NI 

VALU^ 

• 

• 

value  of  t  or  yxv  at  the  end  of  the 
given  history  segment 

41 

ltyp5 

= 

0  ♦  x  is  specified 

1  -*■  ^  is  specified 

42-50 

VALUj 

• 

• 

value  of  x  or  y  at  the  end  of  the 
given  history  segment 

51  LTYP.  =  0  -*■  t  is  specified 

1  ♦  yJL  is  specified 
yz 

52  -  60  VALU,  s  value  of  t  or  v  at  the  end  of  the 

given  history  segment 

66  -  70  N1NC  :  number  of  increments  into  which  the 

given  history  segment  is  to  be  subdivided 

71-80  SRAT  :  incrementing  ratio  which  defines  the 

relative  magnitude  of  two  consecutive 
loading  increments  within  the  given 
history  segment.  A  value  of  SRAT  =  1.0 
results  in  N1NC  equally  spaced  loading 
increments.  By  definition  SRAT  = 
Aon/Aon_I(or  Aen/Aeft-1),  where  the 
superscripts  n-1  and  n  denote  the  values 
in  two  consecutive  loading  increments 
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Experimental  Data  and  Plotting  Instructions  Cards: 


The  following  sequence  of  cards  must  be  supplied  for  each  of  the 
NTST  distinct  experimental  tests  for  which  model  predictions  are 
required: 


1st  Card  (213. 

E10.3,  15): 

Columns 

1  -  5 

ITST 

• 

• 

test  number 

6-10 

NPLT 

• 

• 

number  of  distinct  experimental  relations 
associated  with  the  given  test 

11-20 

WTST 

• 

• 

weighting  factor  for  the  given  test.  By 
default,  WTST  =  1.0 

21-25 

NWPT 

• 

• 

number  of  distinct  experimental  obser¬ 
vations  associated  with  the  given  test  to 
which  weights  other  than  1.0  will  be 
assigned 

2nd  Card  (*  (215.  E10.3)): 

Columns 

1  -  5 

KPLTj 

• 

• 

type  of  plot  for  the  first  experimental 
relation  (see  Table  2.2) 

o 

-H 

1 

NO 

nexpa 

• 

• 

number  of  distinct  experimental  obser¬ 
vations  associated  with  the  given  relation 

11-20 

WPLTj 

• 

• 

weighting  factor  for  the  given  relation. 
By  default,  WPLT  =  1.0 

21  -  SO 

♦  ♦ 

• 

• 

repeat  entries  for  columns  1-20  as 
described  above  (in  the  identical  format) 
until  all  NPLT  distinct  experimental 
relations  have  been  described 

The  following  sequence  of  cards  must  be  provided  for  each  of  the 
NPLT  distinct  experimental  relations  associated  with  test  ITST  for 
which  NEXP  k  0: 

3rd  Cards  (SE10.3): 

Use  as  many  cards  as  necessary  to  enter  the  y  -  values  (or  ordinates) 
of  the  NEXP  distinct  experimental  observations  associated  with 
relation  KPLT  of  test  ITST 

♦th  Cards  (SE10.3): 

Use  as  many  cards  as  necessary  to  enter  the  x-values  (or  abscissae) 
of  the  NEXP  distinct  experimental  observations  associated  with 
relation  KPLT  of  test  ITST 


Special  Experimental  Weightings  Cards:  (include  only  if  £  NWPT.  £ 
0  in  section  VI):  i=l  1 

1st  Cards  (2D.  1».  E10.3): 


Columns 


test  number 


7-10 


11  -  20  WPNT 


21-80  ♦  ♦ 


plot  number 

experimental  observation  number 

weighting^  factor  for  the  observation 

in  the  IP1"  plot  of  the  ITtn  test 

repeat  entries  for  columns  1-20  as 
described  above  (in  the  identical  format) 
until  all  £  NWPT  distinct  special 
weightings  have  been  described. 
Continue  on  a  2nd  card  if  necessary 


The  remaining  input  cards  are  only  necessary  for  those  analyses  where 
JOPT  s  1  (in  section  0)  and,  therefore,  model  calibration  is  required. 


Calibration  Control  Codes  Card: 

1st  Card  (913b 
Columns 

1-5  NDH4  :  number  of  parameters  whose  optimal 

values  are  being  sought 

6-10  NFMX  :  maximum  allowable  number  of  objective 

function  evaluations 

11-15  IPRN  =  1  -*■  only  final  results  printed 

2  *  results  printed  after  each  line 

search 

3  ♦  results  printed  after  each  objective 

function  evaluation 

16  -  20  ICON  =  1  ♦  convergence  based  on  parameter 

values 

2  -*•  convergence  based  on  function 

values 

3  -*•  convergence  based  on  gradient 

values 

21  -  25  1CUB  =  0  ♦  strict  line  search  convergence 

criterion  used 

1  ♦  normal  line  search  convergence 
criterion  used 

26  -  30  IHES  =  0  ♦  Hessian  conditioning  not  performed 

1  -►  Hessian  conditioning  performed 

31  -  35  NHE5  :  number  of  non-zero  values  to  be  entered 

for  the  initial  inverse  Hessian  matrix.  If 
NHES  =  0,  the  inverse  Hessian  will 
initially  be  set  to  the  identity  matrix 

36-40  ISCA  =  0  ♦  parameter  scaling  not  performed 

1  ♦  parameter  scaling  performed 

41  -  45  NSLN  =  0  -*■  no  initial  estimates  computed 

1  ♦  initial  estimates  computed  for 

either  R  or  R 

2  *  initial  estimates  computed  for 

both  R  and  R 


I  -  10  EPFA  :  convergence  criterion  associated  with  the 

objective  function,  f.  If  ICON  =  2  (see 
section  VIII),  the  search  will  be  ter¬ 
minated  when  the  k-th  line  search  has 
failed  to  reduce  f  by  an  appreciable 
amount: 


EPFA  >  «k_,-V 

11-20  EPSA  :  convergence  criterion  associated  with  the 

gradient  of  the  objective  function,  g.  If 
ICON  =  3  (see  section  VIII),  the  search 
will  be  terminated  when  the  slope  of  the 
k-th  computed  search  direction  is  suit¬ 
ably  small: 

EPSA  >  (-gjHkg|<) 

21  -  30  Cl  :  value  of  the  line  search  exclusion  para¬ 

meter  o  A  value  of  a*  0.9  is  recom¬ 
mended 

31  -  40  C2  :  value  of  the  line  search  exclusion  para¬ 

meter  p.  A  value  in  the  range  0.0001 
<  p  <  0.01  is  suggested 

41-30  C3  :  value  of  the  bracket  check  parameter 

t.  A  value  of  t  =  0.1  is  recommended 

51-60  EPXL  :  line  search  exit  criterion  associated  with 

the  parameters  x.  A  value  in  the  range 
0.001  <  EPXL  <  0.01  is  suggested 

61  -  70  EPFL  :  global  fail-safe  exit  criterion  associated 

with  the  objective  function  f 

71-75  NFLX  :  maximum  allowable  number  of  objective 

function  evaluations  that  may  be  made 
before  the  fail-safe  termination  criterion 
(based  on  EPFL)  is  invoked 
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X.  Finite-Difference  Specifications  Card; 


1st  Card  (5E10.3); 
Columns 


1  -  10 

ETAF 

• 

• 

relative  error  associated  with  the 
computation  of  objective  function  values 

11  -  20 

ETAX 

• 

• 

relative  error  associated  with  the  re¬ 
presentation  of  parameter  values  in  a 
finite  word-length  machine 

21-30 

DMIN 

• 

• 

smallest  differencing  interval  permitted 
in  calculating  first  derivatives  by  finite 
differences 

31  -  40 

FDER 

• 

maximum  permissible  relative  truncation 
error.  If  the  estimated  error  is  greater 
than  FDER,  central  rather  than  forward 
differences  are  used  to  compute  the 
gradient.  A  value  in  the  range  0.0001 
<  FDER  <  0.01  is  suggested 

41-50 

SFUN 

• 

• 

scaling  factor  for  the  objective  function. 

The  optimization  code  seeks  the 
minimum  of  SFUN»f.  By  default,  SFUN 
=  1.0 


XL  Initial  Estimates  Card: 

1st  Cards  (8E1Q.3); 

Use  as  many  cards  as  necessary  to  enter  initial  estimates  for  the 
NDIM  parameters  whose  optimal  values  are  being  sought.  For 
typical  values,  see  Table  2.3 


XII.  Lower  Bounds  Card: 

1st  Cards  (8E10.3): 

Use  as  many  cards  as  necessary  to  specify  the  lower  bounds  fa.  < 
x.)  associated  with  each  of  the  NDIM  parameters  whose  optimal 
values  are  being  sought.  For  typical  values,  see  Table  2.3 


XID.  Upper  Bowds  Card; 

1st  Cards  (8E10.3): 

Use  as  many  cards  as  necessary  to  specify  the  upper  bounds  (x.  < 
u.)  associated  with  each  of  the  NDIM  parameters  whose  optimal 
values  are  being  sought.  For  typical  values,  see  Table  2.3 


XIV.  Initial  Differencing  Intervals  Card: 

1st  Cards  PE1Q.3 h 

Use  as  many  cards  as  necessary  to  specify  the  initial  finite 
difference  differencing  intervals  associated  with  each  of  the  NDIM 
parameters  whose  optimal  values  are  being  sought.  For  typical 
values,  see  Table  2.3 


XV.  Scaling  Factors  Cardb  (include  only  if  ISCA  =  1  in  section  Vm):* 
1st  Card  (SE10.3): 

Use  as  many  cards  as  necessary  to  specify  the  scaling  factors 
associated  with  each  of  the  NDIM  parameters  whose  optimal  values 
are  being  sought 


XVL  Convergence  Specifications  Card:  (include  only  if  ICON  =  1  in  section 
VDI): 

1st  Card  (8E10.3): 

Use  as  many  cards  as  necessary  to  specify  the  convergence  criteria 
associated  with  the  parameters  x.  If  ICON  =  1,  the  search  will 
be  terminated  when  the  kth  line  search  has  established  the  location 
of  the  minimum  to  a  sufficiently  high  degree  of  accuracy: 

>  |x.k_1  -  xk|  0=1,  NDIM) 
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XVIL  Initial  Inverse  Hessian  Specification  Card:  (include  only  if  NHES  >  0  in 
section  VIII): 


1st  Card  (1  (213 

.  Eia3))s 

Columns 

1  -  5 

IR 

• 

• 

row  number 

6-10 

IC 

• 

• 

column  number 

11  -  20 

H2LJ 

• 

• 

value  of  the  (IR,  IC)'*1  component  of 
initial  inverse  Hessian  matrix 

21  -  80 

♦  4 

• 

• 

repeat  entries  for  columns  1-20  as 
described  above  On  the  identical  format) 
until  all  NHES  components  have  been 
specified.  Continue  on  a  2nd  card  if 
necessary 

Analytic  Initial 

Estimates  Card: 

(include  only  if  NSLN  >  0  in  section 

vnh 

The  following  card  must  be  repeated  NSLN  times  —  once  for  each 
initial  estimate  required: 

1st  Card  (5E10.3.  I3h 

Columns 

1  -  10 

X 

• 

• 

slope  of  the  isotropic  consolidation  line 
in  e  -  In  p  space 

11  -  20 

IC 

• 

• 

slope  of  the  swell/recompression  line  in 
e  -  In  p  space 

21-30 

TOON 

• 

• 

initial  value  of  the^  mean  normal 
effective  stress  p"  =  p' 

V*) 

1 

o 

PFAL 

• 

• 

critical  value  of  the  mean  normal 
effective  stress  p"  =  p' 

41-50 

RI 

: 

initial  estimate  for  the  surface  shape 
parameter  R.  By  default,  RI  =  2.50 

51  -  55 

KR 

: 

1  ♦  R  *  R  is  assumed 

2  ■*>  R  =  R^  is  assumed 

For  most  analyses,  ISC  A  =  0  and  NHES  =  0,  and  therefore  input  cards 
XV  and  XVII  will  normally  be  omitted. 


'H  1,-TI,  0 


h"lo/R(a)  T_6 


Figure  2.1: 


l0(») 


A  Schematic  Illustration  Of  The  Bounding 
Surface  and  Radial  Mapping  Rule  In 
Invariant  Stress  Space 


TABLE  2.1 


! 

I 


I 


» 


Correspondence  Between  Parameter  Type  and  Property  Number 


Property 

Number 

Parameter 

Type 

Property 

Number 

Parameter 

Type 

1 

X 

11 

T 

2 

< 

12 

VRc 

3 

Mc 

13 

v\ 

4 

Me/Mc 

14 

c 

5 

G 

15 

s 

6 

r 

16 

hc 

7 

Pi 

17 

mc 

8 

Pa 

18 

he/hc 

9 

Rc 

19 

metac 

10 

Ac 

TABLE  2.2 

Correspondence  Between  Relation  Type  and  Relation  Number 

Relation 

Relation 

Number 

Type 

1 

q  vs  p' 

2 

q  vs  Ej 

3 

P'  VS  Ej 

4 

U  VS  E  | 

5 

Ev  VS  €1 

TABLE  2.3 


5 

i 


i 

i 
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Typical  Parameter  Values,  Bounds  and  Initial  Differencing  Intervals 


Parameter 

Number 

Parameter 

Typical 

Value 

Typical 

Lower 

Bound 

Typical 

Upper 

Bound 

Typical  Initial 
Differencing 
Interval 

1 

X 

0.20 

0.10 

0.40 

0.02 

2 

K 

0.04 

0.02 

0.0S 

0.004 

3 

Mc 

1.00 

0.75 

1.25 

0.10 

4 

M  /M 
e  c 

1.00 

0.75 

1.25 

0.05 

5 

G 

2000 

1000 

10000 

200 

7 

Pj. 

pa 

°*25*Pa 

pa 

-0.25«pa 

9 

Rc 

2.50 

2.00 

3.00 

0.05 

10 

Ac 

0.10 

0.03 

0.20 

0.01 

11 

T 

0.10 

0.05 

0.15 

0.05 

12 

VRc 

1.00 

0.75 

1.25 

0.02 

13 

VAc 

1.00 

0.50 

2.00 

0.10 

14 

c 

0.50 

0.00 

0.75 

0.05 

15 

s 

1.00 

1.00 

2.00 

0.50 

16 

h 

c 

0.25 

0.05 

2.00 

0.05 

17 

mc 

0.50 

0.20 

1.00 

0.50 

IS 

he/hc 

1.00 

0.50 

4.00 

0.20 

19 

me/mc 

1.00 

0.50 

4.00 

1.00 

m.  EXPLANATORY  NOTES  AND  INPUT  RECOMMENDATIONS 
LI)  Input  Sections  I  and  II 

The  program  MODCAL  may  be  used  to  accomplish  model  calibration 
,30PT=1),  or  to  generate  model  predictions  (3RUN=1),  or  both.  If  the  code  is 
used  to  generate  model  predictions,  plots  of  the  experimentally  observed  and 
theoretically  computed  response  relations  may  be  requested  (3PLT=1)  to 
sifiplement  the  numerical  output  data. 

The  constants  KRSL  and  KNRM  partially  define  the  manner  by  which  the 
objective  function  is  formed  (see  section  1.2).  If  model  calibration  is  not 
required,  the  values  of  KRSL,  KNRM,  NOPT  and  KOPT(II)  may  be  omitted.  The 
absolute-Euclidean  measure  of  error  (KRSL=0  and  KNRM=0)  has  been  found  to 
produce  the  most  easily  minimized  objective  function.  Hence,  in  most  applications 
this  measure  should  be  employed.  The  vertical  measure  of  error  (KRSL=1)  should 
only  be  used  if  all  of  the  response  relations  included  in  the  calibration  data 
base  (see  section  3.4)  are  either  monotonically  increasing  or  monotonically 
decreasing  with  respect  to  their  x-values  or  abscissae.  For  example,  the  vertical 
measure  could  be  used  with  curves  a-g  of  Figure  3.1,  but  not  with  curve  h 
(since  this  relation  first  increases  and  then  decreases  with  respect  to  its 
x-coordinates).  Thus,  the  vertical  measure  should  normally  not  be  used  if  q  vs 
p'  relations  are  included  in  the  calibration  data  base. 

The  NOPT  parameters  to  be  established  by  optimization  are  identified  by 
entering  their  respective  "property  numbers"  in  the  KOPT(H)  array.  The 
correspondence  between  parameter  type  and  property  number  is  defined  in  Table 
2.1.  As  may  be  seen,  the  ordering  scheme  used  in  Table  2.1  is  identical  to 
that  employed  in  input  section  III;  that  is,  X  is  property  1,  <  is  property  2,  and 
so  forth. 
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3.2)  Input  Section  in 

The  values  of  the  required  material  and  model  properties  are  specified 
in  input  section  ID.  A  more  detailed  description  of  the  various  parameters, 
including  an  examination  of  their  qualitative  and  quantitative  influence  on  the 
predicted  material  response,  is  presented  by  DeNatale  (1982).  If  model  calibration 
is  required,  it  is  not  necessary  to  enter  estimates  for  those  parameters  whose 
optimal  values  are  being  sought. 

3.3)  Input  Sections  IV  and  V 

As  mentioned  earlier  in  section  1.4,  the  code  relies  on  two  subroutines 
— EVAL  and  CLAY  —  to  generate  the  required  model  predictions.  These 
predictions  can  be  acquired  by  means  of  either  a  "conventional"  (KIND=  1)  or  a 
"reformulated'  (KIND=0)  analysis.  The  reformulated  analysis  is  a  modification 
of  the  mixed  finite  element  procedure  developed  by  Herrmann  (1965)  for 
incompressible  and  nearly  incompressible  solids.  In  it,  the  six  strain  components 
together  with  the  pore  water  pressure  represent  the  "mixed"  set  of  primary 
dependent  variables.  When  both  the  conventional  and  reformulated  analyses 
converge,  they  yield  identical  results.  However,  for  certain  undrained  and 
near-failure  conditions  the  reformulated  analysis  may  successfully  converge,  while 
the  conventional  analysis  may  not.  Additional  details  of  the  two  formulations 
are  provided  by  Herrmann  et  al  (1981,  1982). 

Both  types  of  analyses  use  the  method  of  successive  approximation,  and 
acceleration  factors  can  be  applied  to  the  stress  and  strain  vectors  to  increase 
the  rate  of  convergence.  The  acceleration  factors  are  calculated  by  means  of 
the  Aitken's  formula  (see,  for  example,  Isaacson  and  Keller  (1966)),  and  the 
maximum  and  minimum  factor  values  are  controlled  by  the  constant  CNFR.  By 
default,  CNFR=0.3,  and  if  CNFR=1.0,  acceleration  factors  will  not  be  used. 


Finally,  the  analyses  may  be  accomplished  using  either  engineering 
(LARG=0)  or  natural  (LARG=1)  strains.  This  option  was  introduced  solely  to 
facilitate  the  comparison  of  model  predictions  with  experimental  observations 
reported  in  either  format. 

In  order  to  simulate  a  given  experimental  test  ITST,  the  applied  loading 
history  must  be  specified  in  piecewise  form  by  entering  descriptions  of  its  NSEG 
distinct  "history  segments."  A  particular  history  segment  is  defined  by  the 
values  of  the  six  stress  and/or  strain  components  (LTYPi  and  VALlh;  i=l,6) 
which  act  on  the  material  element  at  the  end  of  that  segment.  Within  the 
program,  the  segment  is  further  subdivided  into  NINC  separate  increments,  and 
iteration  is  performed  within  each  increment.  The  relative  size  of  the  increments 
which  make  up  a  given  history  segment  is  controlled  by  the  constant  SRAT.  By 
default  SRAT=1.0,  and  the  segment  is  then  subdivided  into  NINC  equal  sized 
increments.  If,  in  any  particular  increment,  convergence  (as  defined  by  the 
constant  ERMX)  is  not  achieved  within  ITMX  iterations,  an  appropriate  message 
is  printed  and  the  analysis  will  then  stop.  By  default,  ERMXrO.Ol  and  ITMX=10. 

A  conventional  monotonic  triaxial  loading  could  possibly  be  described  with 
as  few  as  one  history  segment.  For  example,  a  strain-controlled  test  which 
began  at  an  all-around  confining  stress  of  o^.  =  oq  and  ended,  after  30  equal 
sized  steps,  at  an  axial  strain  of  Cj  =  15%,  could  be  simulated  by  specifying: 

°x  =  °o 

O  =0 

y  o 

e2  =  0.15 

v  =  Y  =  Y  =  0.0 
•xy  'xz  ’yz 

NINC  =  30 
SRAT  =  1.0 
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There  are  three  situations  where  it  would  be  necessary  to  use  more  than  one 
segment  to  describe  a  particular  loading  history;  namely: 

i)  if  one  or  more  of  the  six  loading  components  changes  during  the 
test  from  a  stress  quantity  to  a  strain  quantity  (or  vice  versa); 

ii)  if  the  loading  history  cannot  be  subdivided  into  increments  of  the 
desired  length  through  the  use  of  the  parameters  NINC  and  SRAT, 
alone;  or, 

iii)  if  a  cyclic  loading  history  must  be  described. 

For  example,  a  minimum  of  two  segments  would  be  required  if  the  triaxial 
specimen  previously  described  was  first  loaded  to  an  axial  strain  of  e  =15%  and 
then  unloaded  to  its  initial  stress  state  Oy.  =  oq  in  10  equal  sized  steps: 


e,  *  0.15 

Z 


o  =  o 

y  o 

cr,  =  a 

z  o 


Yxy  ~  ^xz  “  ^yz 
NINC  =  30 


Yxy  =  Yxz  Yyz 
NINC  =  10 


SRAT  =  1.0 


SRAT  =  1.0 


Or,  if  only  loading  was  involved,  but  it  was  desired  to  use  increments  of  Aez 

=  0.2%  for  axial  strains  of  e_  <  2%  and  increments  of  Ae  =  1.0%  for  axial 

z  z 

strains  of  e  >  2%,  a  minimum  of  two  segments  would  again  be  required: 

Z 


°x  =°o 


O  =0 

y  o 

e  =  0.02 

Z 


v  s  y  =  Y  =  0.0 
Txy  ’xz  lyz 

NINC  =  10 


(ii)  ox.o0 
°y  =  °o 


e,  =  0.15 

Z 

V  =  Y«  =  V 

NINC  =  13 


SRAT  =  1.0 


SRAT  =  1.0 
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Careful  thought  should  be  given  to  the  total  number  of  increments  used 
to  represent  a  particular  loading  history.  The  importance  of  increment  size  in 
EVAL-directed  analyses  is  analogous  to  the  significance  of  grid  size  in  a  general 
finite  element  analysis;  that  is,  the  predicted  material  response  may,  to  some 
extent,  be  influenced  by  the  size  of  the  increments  into  which  the  applied 
loading  history  is  subdivided.  If  too  few  increments  are  used,  convergence  may 
not  occur,  or  the  program  may  converge  to  an  "incorrect"  state  of  stress  and/or 
strain.  If  too  many  increments  are  used,  the  analysis  would  become  unnecessarily 
expensive.  Note  that  the  cost  of  a  given  analysis  is  controlled  not  by  the 
number  of  history  segments  or  loading  increments  used,  but,  rather,  by  the  total 
number  of  iterations  required,  ail  increments  considered.  Thus,  an  analysis  which 
uses  a  greater  number  of  increments  may  actually  be  more  economical,  if  it 
can  be  completed  in  a  lesser  total  number  of  iterations. 

Typically,  from  20-40  increments  are  sufficient  to  describe  a  particular 
monotonic  loading.  In  general,  undrained  analyses  require  smaller  increments 
than  drained  analyses,  and  the  unloading  portions  of  a  test  can  normally  be 
described  with  larger  increments  than  the  loading  segments.  The  code  is  designed 
to  simulate  tests  which  begin  at  a  hydrostatic  state  of  stress  and 

relatively  small  increments  should  be  used  if  and  when  the  stress  state  first 
moves  off  the  hydrostatic  (Ij)  axis.  An  automated  increment  generation  scheme 
is  currently  under  development  as  part  of  the  ongoing  research  with  the  Bounding 
Surface  model.  However,  since  the  scheme  has  not  yet  been  incorporated  into 
the  computer  code,  some  user-initiated  checks  should  first  be  made  with  respect 
to  optimal  increment  size,  before  proceeding  with  any  potentially  expensive 
analyses. 

The  parameters  PLIM  and  TLIM  are  the  second  group  of  constants 
associated  with  the  formation  of  the  objective  function.  Their  roles  can  perhaps 
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best  be  understood  by  considering  the  expression  for  residuals  used  in  the  code. 
The  unweighted  residual  corresponding  to  experimental  observation  k  is  given 


where: 


rk  yD 

•  5? 


yu 

If  (~  >  PLIM), 
YP 

If  (-*  <  PLIM), 
yP 


yl  =  yk 


Yj  =  yp  •  PLIM 


If  (-2  >  TLIM), 
yr 

yn 

If  (-E  <  TLIM), 
yr 

and  where  (refer  also  to  Figure  3.2): 


y2  =  yp 

y2  =  yr  •  tlim 


rk  =  the  normalized  distance  between  the  kth  experimental 
observation  and  the  prediction  curve; 
yk  =  the  y-value  (or  ordinate)  of  the  kth  experimental  obser¬ 
vation; 

yp  =  the  maximum  absolute  experimentally  observed  y-value, 
all  points  of  the  given  relation  consider  ?d; 
yr  =  the  maximum  absolute  experimentally  observed  y-value, 
all  relations  of  the  given  type  considered; 

Yj  =  the  base  y-value  for  residual  scaling  for  all  points  of 
a  given  relation;  and, 

Yj  =  the  base  y-value  for  residual  scaling  for  all  relations 
of  a  given  type. 


As  pointed  out  earlier  in  section  1.2,  before  an  objective  function  can  be  formed, 
the  user  must  define  to  the  code  what  is  meant  by  "equal"  errors.  If  equality 
is  based  on  relative  errors  throughout,  then: 


PLIM  =  0 
TLIM  =  0 


*1  =  yk 

y2  =  yp 


and  therefore  the  residual  R^  will  be  the  same  at  two  points  A  and  B  only  if 
there  are  identical  relative  errors  (r^/y^)  at  A  and  B.  If  equality  is  based  on 
absolute  errors  throughout,  then: 

PLIM  r 
TLIM  = 

yl  = 
y2  = 


and  therefore  the  residual  will  be  the  same  at  two  points  A  and  B  only  if 
there  are  identical  absolute  errors  r^  at  A  and  B.  Additionally,  equality  could 
be  based  on  relative  errors  within  a  given  relation  (PLIM=0)  but  absolute  errors 
between  relations  of  a  given  kind  (TLIM=1),  or  vice  versa  (PLIM=1  and  TLIM=0), 
or  anywhere  in  between  (0  <  PLIM  <  1  and  0  <  TLIM  <  1).  By  default, 
PLIM=0.01  and  TLIM=0.01. 


The  most  serious  argument  against  the  use  of  relative  errors  is  that  small 
deviations  at  low  magnitude  observations  can  produoe  extremely  large  contribu¬ 
tions  to  the  total  computed  overall  error.  In  other  words,  as  yk  approaches 
zero,  the  quantity  (r^/y^)  approaches  infinity,  regardless  of  the  absolute  error 
r^.  The  constants  PLIM  and  TL1M  permit  the  use  of  relative  errors  over  most 
of  the  data  base,  while  preventing  those  observations  whose  magnitudes  are 
relatively  small  from  wduly  influencing  the  computed  overall  "goodness"  of  the 
model  predictions.  The  net  effect  of  a  given  combination  of  PLIM  and  TLIM 
could  probably  be  duplicated  through  a  judiciously  chosen  set  of  weighting  factors 
in  input  section  VI  (see  section  3.4).  However,  weighting  factors  were  introduced 
for  a  distinctly  different  reason,  as  will  be  considered  in  section  3.4. 

3.4)  Input  Sections  VI  and  VD 

The  type  of  plots  associated  with  a  given  test  ITST  are  identified  by 
entering  their  respective  "relation  numbers"  as  KPLT.  The  correspondence 
between  plot  type  and  relation  number  is  defined  in  Table  2.2.  As  may  be  seen, 
a  q  vs  p'  plot  is  assigned  a  relation  number  of  KPLT=1,  and  so  forth. 

If  NEXP/O,  the  individual  data  points  associated  with  the  experimental 
response  relation  KPLT  must  be  entered  in  the  same  order  in  which  they  were 
originally  observed,  and  not  in  some  random  fashion.  In  other  words,  the  first 
point  specified  must  correspond  to  the  initial  experimental  measurement,  while 
the  last  point  entered  must  represent  the  final  observation.  This  ordering  scheme 
must  hold  for  both  the  3rd  and  4th  card  group  entries  of  input  section  VI.  If 
NEXP=0,  the  3rd  and  4th  cards  associated  with  the  particular  relation  KPLT 
must  be  omitted. 

The  weighting  factors  WTST,  WPLT  and  WPNT  are  the  third  and  final 
group  of  constants  associated  with  the  formation  of  the  objective  function.  Like 
KRSL  and  KNRM  (see  section  3.1)  and  PLIM  and  TLIM  (see  section  3.3),  their 


values  may  be  omitted  if  model  calibration  is  not  required.  By  default,  all 
tests,  relations  and  individual  observations  are  weighted  equally.  However, 
different  weights  may  be  assigned  to  specific  components  of  the  data  base  if 
it  is  felt  that  certain  tests,  relations  or  observations  are  more  reliable  or 
representative  than  others,  or  if  it  is  necessary  to  have  the  final  model  predictions 
fit  some  data  more  closely  than  others.  Thus,  by  assigning  a  weight  of  WPNT=0 
to  a  given  experimental  observation,  it  becomes  possible  to  include  that  point 
in  the  calibration  data  base,  without  it  influencing  the  computed  value  of  the 
objective  function. 

If  model  calibration  is  required,  the  content  of  the  calibration  data  base 
should  be  complete  and  diverse  enough  to  permit  each  of  the  unknown  parameters 
to  be  uniquely  determined.  The  recommended  calibration  data  base  is  described 
in  detail  in  section  1.5,  and  the  consequences  of  an  inadequate  or  incomplete 
data  base  are  considered  in  section  1.6. 

3.5)  Input  Sections  Vm  and  DC 

As  mentioned  previously  in  section  1.3,  the  line  search  strategy  is  an 
integral  part  of  the  global  optimization  routine.  The  basic  structure  of  the  line 
search  algorithm  is  illustrated  in  Figure  3.3.  As  may  be  seen,  when  the  strict 
line  search  termination  criterion  is  employed  (ICUB=0),  an  additional  situation 
is  introduced  wherein  cubic  interpolation  is  used  to  reduce  the  size  of  the  search 
bracket.  The  primary  advantage  of  this  strict  termination  criterion  is  that  it 
permits  exact  line  searches  to  be  carried  out  in  the  limit  as  o  approaches  zero. 
However,  in  most  practical  problems  inexact  line  searches  (say  a  *  0.7-0.9) 
normally  result  in  a  lesser  total  number  of  objective  function  evaluations.  Hence, 
provided  that  a  a  greater  than  about  0.7,  ICUB=0  and  ICUB=1  analyses  will 
generally  yield  equally  efficient  results. 


The  constants  Cl,  C2  and  C3  represent  the  line  search  parameters  o,  p 
and  x,  respectively.  As  indicated  in  Figure  3.3,  the  parameters  a  and  p  control 
the  degree  of  accuracy  to  which  a  particular  line  search  is  carried  out;  observe 
that  as  a  increases  and  p  decreases,  the  line  search  becomes  less  and  less  exact. 
The  parameter  t  prevents  evaluations  from  being  made  at  points  which  are 
trivially  close  to  the  boundaries  of  the  current  search  bracket.  By  default,  a 
=  0.90,  p  =  0.0001  and  T  =  0.10. 

Normally,  the  line  search  terminates  because  the  exit  criterion  shown  in 
Figure  3.3  has  been  satisfied.  However,  near  the  solution  (that  is,  the  minimum) 
round-off  errors  may  become  significant,  and  it  is  sometimes  more  efficient  and 
reliable  to  use  parameter  values  as  the  exit  criterion.  Thus,  in  the  code,  the 
line  search  is  terminated  whenever  the  predicted  location  of  the  optimum  is  no 
longer  changing  by  a  significant  amount.  Mathematically,  this  additional  exit 
criterion  is  given  by: 


EPXL 


> 


(i  =  1,  ISDIM) 


where  x‘  j  and  x1  represent  the  values  of  parameter  x1  at  two  consecutive 
objective  function  evaluations.  By  default,  EPXL=0.001. 

There  are  a  variety  of  criteria  that  may  be  used  to  test  for  global 
convergence  and  which  may  therefore  serve  as  a  basis  for  terminating  the 
optimization  search.  Three  of  the  most  commonly  used  and  effective  criteria 
are: 


i)  based  on  parameter  values  (!CON=l), 


ii)  based  on  function  values  (ICON =2), 

EPFA  >  -  fk> 

iii)  based  on  gradient  values  (lCON=3), 

EPSA  >  l-g^H^I 

where  the  subscripts  k-1  and  k  refer  to  the  value  of  the  quantity  at  the 
completion  of  the  k-l^  and  kth  line  search,  respectively. 

The  use  of  criterion  (i)  ensures  that  the  location  of  the  optimum  has  been 
determined  to  within  the  specified  accuracy,  but  does  not  guarantee  that  the 
smallest  function  value  has  been  found.  The  use  of  criterion  (ii)  requires  the 
user  to  identify  on  input  the  maximum  absolute  error  that  can  be  tolerated  with 
respect  to  the  magnitude  of  the  objective  function  at  the  optimum.  The  use 
of  criterion  (iii)  may  perhaps  provide  the  best  results;  when  the  gradient  is  small 
in  the  next  computed  search  direction,  any  additional  searches  in  that  direction 
should  yield  little  or  no  reduction  in  the  value  of  the  objective  function. 

Regardless  of  which  of  the  above  three  global  termination  criteria  is  used, 
a  fail-safe  exit  condition  is  also  incorporated,  to  prevent  the  algorithm  from 
continuing  when  no  substantial  progress  is  being  made.  If  at  some  point  in  the 


search 


EPFL  >  <fk.,-ik> 


and  this  condition  continues  to  be  satisfied  for  the  next  NFLX  objective  function 
evaluations,  the  search  is  terminated,  regardless  of  ICON.  This  fail-safe  exit 
condition  was  introduced  to  prevent  unnecessary  function  evaluations  from  being 
made  when  either  parameter  values  (ICON=l)  or  gradient  values  (ICON=3)  are 
used  as  the  convergence  criterion,  but  where  this  criterion  is,  for  one  reason 
or  another,  too  strict  to  ever  be  invoked.  If  function  values  are  used  as  the 
convergence  criterion  (ICON=2),  the  user  should  set  EPFA=EPFL. 


The  codes  IHES  and  ISCA  (and  the  constant  SFUN  of  input  section  X) 
are  used  to  select  options  which  may  improve  the  conditioning  of  the  problem 
and  therefore  the  efficiency  of  the  search.  If  IHES=0,  the  initial  inverse  Hessian 
will  be  scaled  prior  to  the  first  update.  This  scaling  procedure  was  originally 
suggested  by  Shanno  and  Phua  (197S)  as  a  means  for  reducing  the  cancellation 
errors  which  may  occur  when  the  terms  in  the  updating  formula  (see  section 
1.3)  are  of  significantly  different  magnitudes.  However,  the  results  of  the 
procedure  have,  in  practice,  been  rather  mixed  (see  Shanno  and  Phua  (1978)). 
For  small  to  medium  sized  problems  (say  NDIM  <  10  to  15),  such  as  those 
normally  encountered  in  model  calibration,  it  is  probably  best  to  set  IHES=0, 
and  thereby  prevent  scaling.  For  larger  problems,  it  may  be  possible  to  improve 
the  algorithm's  performance  by  setting  IHES=1. 

When  the  typical  values  of  the  constitutive  parameters  are  of  vastly 
different  magnitudes  (such  as  a  shear  modulus  of  G=2500  psi,  versus  a  projection 
point  of  c=0.25),  the  Quasi-Newton  strategy  may  perform  quite  poorly  as  a  result 
of  badly  scaled  x,  g  and  6  vectors  (see  section  1.3).  It  is  often  possible,  though, 
to  improve  the  routine's  performance  by  making  the  variable  transformation 

y  =  Dx 

and  minimize,  instead,  f(x(y)]  =  f(y),  where  D  is  a  diagonal  matrix  consisting 
of  NDIM  scaling  factors.  However,  since  the  program  already  employs  a 
trigonometric  variable  transformation  to  enable  the  imposition  of  upper  and 
lower  bounds  (see  section  1.3),  which  generally  results  in  a  very  well  scaled 
problem,  no  additional  scaling  should  normally  be  required  (ISCA=0). 

Finally,  it  is  sometimes  possible  to  improve  the  algorithm's  performance 
by  minimizing  a  scalar  multiple  of  the  objective  function  (SFUN  •  f),  rather 
than  the  function  itself.  This  phenomenon  again  appears  to  be  due  entirely  to 
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the  presence  of  round-off  errors.  However,  since  the  optimal  value  of  the  scalar 
SFUN  is  extremely  dependent  on  the  specific  problem  being  considered,  there 
is  normally  no  reason  to  initially  invoke  the  function  scaling  option  (and  thus, 
by  default,  SFUN=1.0). 

3.6)  Input  Section  X 

As  mentioned  previously  in  section  1.3,  gradients  are  evaluated  by  means 

of  finite  differences,  in  accordance  with  the  strategy  originally  suggested  by 

Stewart  (1967).  To  enable  this  strategy  to  be  efficiently  carried  out,  certain 

problem  and  machine  dependent  tolerance  information  must  be  supplied. 

On  a  finite  word  length  machine,  a  given  property  value  x.  cannot  be 

represented  with  infinite  precision,  but,  rather,  must  be  stored  as  some  truncated 
* 

approximation  x.,  where: 

x*  «  (1  ♦  c j )x .  and  Ej  <  ETAX 

Similarly,  instead  of  providing  a  perfectly  accurate  value  of  the  objective  function 

* 

f,  the  algorithm  can  only  yield  an  approximation  f  ,  where: 

f*  =  (1  +  e2)f  and  e2  <  ETAF 

The  user  must  further  specifiy  the  maximum  truncation  error  that  can  be 
tolerated  in  the  forward  difference  representation  of  derivatives,  FDER.  The 
constant  FDER  thus  permits  the  code  to  identify  when  it  is  necessary  to  replace 
the  forward  difference  formula  with  the  more  expensive  central  difference 
representation.  Finally,  in  order  to  prevent  the  computed  differencing  intervals 
from  becoming  too  small  to  yield  a  representative  gradient,  a  lower  limit  on 
the  size  of  the  intervals  must  also  be  specified  by  means  of  the  constant  DM1N. 
Since  the  interval  size  check  is  made  after  the  trigonometric  variable  transforma- 


tion  discussed  in  section  1.3  has  been  performed,  this  single  bound  DMIN  can 
be  applied  to  all  components  of  the  gradient  vector  and  still  produce  the  desired 
effect.  By  default,  ETAX=10"*,  ETAF=10*4,  FDEFU10"4  and  DMIN=0.005. 

3.7)  Input  Sections  XI  Through  XVI 

The  ordering  of  the  specifications  in  input  sections  XI  through  XVI  should 
be  consistent  with  that  used  in  input  section  II  to  identify  the  unknown  constitutive 
parameters.  For  example,  if  a  four-dimensional  (NDIM=4)  analysis  was  performed 
to  establish  the  optimal  values  of  Rc,  Ac,  c  and  hc  (having  property  numbers 
of  9,  10,  15  and  16,  respectively,  as  shown  in  Table  2.1),  and  if  the  specification 
of  input  section  n  was  made  in  the  order: 

KOPT  (1)  =  9 

KOPT  (2)  =  10 
KOPT  (3)  =  15 
KOPT  <*)  *  16 

then  the  initial  estimate  for  Rc  should  occupy  the  first  position  (columns  1-10) 
of  input  section  XI,  the  lower  bound  for  c  should  occupy  the  third  position 
(columns  21-30)  of  input  section  XII,  and  so  forth.  For  convenience,  typical 
starting  estimates  and  lower  and  upper  bounds  are  listed  in  Table  2.3. 

The  initial  differencing  intervals  of  input  section  XIV  should  be  large 
enough  to  yield  an  appreciable  change  in  the  objective  function,  but  not  so  large 
that  the  local  geometry  of  the  fwction  is  misrepresented.  As  a  rule  of  thumb, 
the  size  of  the  initial  differencing  intervals  should  be  about  5-10%  as  great  as 
either  (i)  the  initial  estimate  for  the  corresponding  parameter,  or  (ii)  the 
difference  between  the  specified  ipper  and  lower  bounds.  For  convenience, 
suggested  interval  sizes  are  given  in  Table  2.3. 


3.8)  Input  Section  XVII 

By  default,  the  initial  inverse  Hessian  matrix  is  set  to  H  =  I  (see  section 
1.3).  However,  if  the  objective  function's  curvature  at  the  starting  location  is 
known,  an  improved  initial  estimate  for  H  could  be  specified  by  setting  NHES 
>  0  in  section  VM  and  entering  the  necessary  second  derivative  information  in 
section  XVD.  Since  the  initial  inverse  Hessian  must  be  both  symmetric  and 
positive-definite,  it  is  only  necessary  to  enter  the  non-zero  components  on  the 
diagonal  and  in  the  upper  (or  lower)  triangular  block. 

3.9)  Input  Section  XVm 

Finally,  if  either  Rc  or  Rg  (or  both)  are  among  the  parameters  whose 
optimal  values  are  being  sought,  improved  initial  estimates  can  be  acquired  by 
setting  NSLN  t  0  in  section  VM  and  entering  the  necessary  data  in  section 
XVID.  These  improved  estimates  are  generated  internally,  by  identifying  the 
roots  of  the  equation  which  defines  the  undrained  stress  path  of  a  normally 
consolidated  soil: 

F(R)  =  0  =  2(oB  -  a28)  +  (a28  -  2)R  4  2R2  -  R3 


where: 


X  -  K 


and  where  X  and  k  represent  the  slopes  of  the  material's  isotropic  consolidation 
and  swell/recompression  curves  in  e  -  In  p"  space,  and  p*  and  p^  denote  the 
initial  and  critical  (or  "failure")  values  of  the  mean  normal  effective  stress  p", 
as  indicated  in  Figure  3.4. 


Mean  Normal  Effective  Stress  p 


Figure  3.#:  The  Undrained  Stress  Path  For  A  Typical 
Consolidated  Soil  Specimen 


Normally 


i)  estimate  the  location  of  the  root  R 


ii)  compute  u(R)  =  f 

iii)  compute  u"(R)  =  1  -  ^ 

[f'(R)r 

iv)  compute  Rn+1  =  Rn  -  =  Rn  +  6° 


v)  if  c  >  |  6n  |  then  stop,  otherwise  repeat  steps  (ii)-(iv) 


Although  a  cubic  equation,  in  general,  possesses  three  roots,  the  equation 
presented  above  has  normally  been  found  to  have  only  a  single  real  root  in  the 
approximate  range  1.5  <  R  <  3.5  (where  the  final  value  of  R  should  nearly 
always  lie).  Hence,  the  internally  computed  estimate  for  R  should  typically 
lead  to  a  rather  good  correlation  between  the  observed  and  predicted  responses 
associated  with  a  normally  consolidated  soil.  If  improved  estimates  for  both 
Rc  and  Rg  are  required  (NSLN  =  2),  the  specifications  pertaining  to  Rc  must 
be  entered  first. 


IV.  MODIFYING  THE  CODE 

The  program  MODCAL  incorporates  29  separate  subroutines,  each  of  which 
performs  a  limited  and  specific  task.  This  modular  form  gives  the  code  maximum 
generality  and  extreme  versatility.  For  example,  if,  in  the  future,  another 
optimization  strategy  is  found  to  be  more  effective  than  the  Quasi-Newton 
procedure  presently  employed,  the  associated  subroutines  may  be  readily  replaced, 
without  having  to  additionally  modify  the  other  components  of  the  code.  Or, 
if  it  was  desired  to  use  boundary-value  measurements  to  calibrate  the  material 
model,  rather  than  homogeneous  test  results,  subroutine  EVAL  could  easily  be 
replaced  with  an  appropriate  multi-element  finite  element  program  (such  as  the 
nonlinear  two-dimensional  code  NTD  developed  by  Herrmann  et  al  (1982b)). 
Finally,  if  it  was  desired  to  adopt  the  present  automated  calibration  strategy 
to  another  constitutive  formulation,  subroutines  CLAY  and  BOUNDS  could  rapidly 
be  replaced  with  the  appropriate  material  model.  The  precise  procedure  for 
linking  a  new  materials  subroutine  to  the  driving  subroutine  EVAL  is  discussed 
in  detail  by  Herrmann  et  al  (1980,  1981). 

The  code's  modular  structure  also  enables  it  to  be  readily  implemented 
on  limited  core  memory  minicomputers,  such  as  the  LSI- 11/23.  The  presence 
of  many  small  subroutines  facilitates  the  use  of  "overlays",  wherein  the  various 
subroutines  are  swapped  into  and  out  of  core  as  they  are  needed.  The  use  of 
overlays  therefore  often  makes  it  possible  to  accommodate  large  programs  which 
would  otherwise  exceed  the  core  memory  capacity  of  a  small  machine. 
Alternatively,  if  the  code  is  to  be  used  only  in  a  limited  form  (such  as  for 
performing  model  calibration  only,  but  without  plotting),  the  unnecessary 
subroutines  may  simply  be  removed,  thus  reducing  the  size  of  the  program. 

In  order  to  facilitate  the  use  of  these  various  options,  and  to  enable  the 
user  to  more  clearly  identify  where  particular  modifications  may  be  required, 
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a  listing  of  the  29  subroutines  and  their  respective  functions  will  now  be 
presented: 


Subroutine 

1.  MODCAL: 

2.  DATAIN: 

3.  FNDMAX: 

4.  EVAL: 

3.  AITKEN: 

6.  PLTCHK: 

7.  JSDPLT: 

5.  LSHIFT: 

9.  BORDER: 

10.  RSDUAL: 


Function 


the  main  driving  program  for  the  Bounding  Surface 
soil  plasticity  model  calibration  and  prediction  code 
directs  the  reading  of  the  required  input  information 
and  control  specifications 

establishes  the  maxima  of  the  specified  experimental 
observations 

directs  the  single-element,  incremental-iterative  finite 
element  analysis  which  generates  the  set  of  model 
predictions  for  the  specified  homogeneous  test  condi¬ 
tions 

computes  the  appropriate  Aitken's  acceleration  factors 
checks  whether  a  given  experimentally  observed  or 
theoretically  predicted  response  relation  is  to  be 
plotted 

directs  the  plotting  of  the  experimental  observations 
and  model  predictions 

adjusts  the  minimum  and  maximum  axes  values,  as 
directed  by  J5DPLT 

produces  the  plot  headings  and  axes  lables,  as  directed 
by  JSDPLT 

directs  the  evaluation  of  the  weighted  residuals  and 
forms  the  governing  objective  function 
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Subroutine 


Function 


11.  EUCLID:  computes  the  minimum  Euclidean  distance  between 

the  experimental  observation  and  the  prediction  curve 

12.  VRT1CL:  computes  the  minimum  vertical  distance  between  the 

experimental  observation  and  the  prediction  curve 

13.  BSMOPT:  directs  the  global  search  component  of  the  model 

calibration  algorithm 

14.  CALDAT:  directs  the  reading  of  all  additional  specifications 

needed  for  model  calibration 

15.  NEWTON:  computes  improved  starting  estimates  for  the  surface 

shape  parameters  Rc  and/or  Rft 

16.  SEARCH:  directs  the  line  search  component  of  the  minimization 

algorithm 

17.  BFGSUP:  computes  the  inverse  Hessian  matrix  H,  in  accordance 

with  the  BFGS  (Broyden-Fletcher-Goldfarb-Shanno) 
updating  formula 

IS.  CONCHK:  checks  whether  the  optimization  algorithm  has 

converged  to  the  minimum 

19.  PRNOUT:  directs  the  printing  of  the  minimization  search  results 

__  *r 

20.  DTPROD:  computes  the  scalar  (dot)  product  of  two  vectors  v  v 

21.  WPROD:  computes  the  product  of  the  vector-vector  operation 

T 

w 

22.  MVPROD:  computes  the  product  of  the  matrix-vector  operation 

Mv 


23.  MMPROD:  computes  the  product  of  the  matrix-matrix  operation 


Subroutine 


Function 


24.  FFUN:  computes  the  value  of  the  objective  function 

25.  GFUN:  computes  the  gradient  of  the  objective  function 

26.  OPEN:  opens  the  required  input  and  output  files 

27.  EXIT:  closes  the  required  input  and  output  files 

28.  CLAY:  contains  the  numerical  implementation  of  the  govern¬ 

ing  constitutive  equations,  and  thus  provides  the  appro¬ 
priate  theoretical  material  response  to  the  given  stress 
and/or  strain  increment,  as  supplied  by  EVAL 
evaluates  the  relationship  of  the  stress  state  to  the 
bounding  surface,  as  required  by  CLAY. 


29.  BOUNDS: 


V.  EXPANDING  THE  CAPABILITIES  OF  THE  CODE 

The  program  MODCAL  employs  a  number  of  1-,  2-  and  3-dimensional 
arrays.  In  order  to  minimize  storage  requirements,  or  to  enable  the  code  to 
accommodate  larger  problems,  it  may  at  times  be  necessary  to  adjust  the 
dimensions  of  the  size-dependent  arrays.  The  list  which  follows  describes  the 
various  size-dependent  arrays,  dimensioned  in  terms  of  a  few  key  constants: 

in  COMMON  BLKA: 

KOPT  (NDIM) 

in  COMMON  BLKBt  j 

¥1  (NTST) 

¥2  (NTST,  NPLT) 

¥3  (NTST,  NPLT,  NEXP) 
in  COMMON  BLKC: 

NSEG  (NTST) 

PRP2  (NTST,  4) 

LTYP  (NTST,  NSEG,  7) 

VALU  (NT5T,  NSEG,  7) 
in  COMMON  BLKD: 

NPLT  (NTST) 

KPLT  (NTST,  NPLT) 

NEXP  (NTST,  NPLT) 

XV  (NTST,  NPLT,  NEXP) 

YV  (NTST,  NPLT,  NEXP) 

PINC  (NINC  +  1,  NVAR) 
in  COMMON  BLKE: 

XMXR  (NRLN) 

YMXR  (NRLN) 


XMXP  (NTST,  NPLT) 
YMXP  (NTST,  NPLT) 
in  COMMON  BLK3: 

XI  (NDIM) 

X2  (NDIM) 
in  COMMON  BLK4; 

XE  (NDIM) 
in  COMMON  BLK5: 

DL  (NDIM) 

D2  (NDIM) 

DG  (NDIM) 

DX  (NDIM) 

G1  (NDIM) 

G2  (NDIM) 

HI  (NDIM,  NDIM) 

H2  (NDIM,  NDIM) 
in  COMMON  BLK7: 

SF  (NDIM) 

XL  (NDIM) 

XU  (NDIM) 
in  MAIN: 

XV  (NDIM) 
in  subroutine  EVAL: 

XOPT  (NDIM) 
in  subroutine  PLTCHK; 
IPLT  (NRLN) 
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in  subroutine  JSDPLT: 


where 


ICEX  (NEXP) 

ICEY  (NEXP) 

ICPX  (NINC  +  1) 

ICPY  (NINC  +  1) 

in  subroutines  BSMOPT,  FFUN  and  GFUN; 

XX  (NDIM) 

in  subroutine  SEARCH: 

DM  (NDIM) 

XM  (NDIM) 

in  subroutine  BFGSUP: 

SV  (NDIM) 

SM  (NDIM,  NDIM) 
in  subroutine  CONCHK: 

XV  (NDIM) 

in  subroutine  PRNOUT; 

XS  (NDIM) 

in  subroutines  DTPROD,  VVPROD.  MVPROD  and  MMPROD: 

VI  (NDIM) 

V2  (NDIM) 

V3  (NDIM,  NDIM) 

V*  (NDIM,  NDIM) 
i  by  definition: 

NDIM  &  number  of  parameters  whose  optimal  values  are  to  be 
established  by  optimization; 

NTST  >  number  of  distinct  experimental  tests; 

NPLT  >  maximum  number  of  response  relations  associated  with  a 


single  experimental  test; 


NEXP  £  maximum  number  of  experimental  observations  associated 
with  a  single  response  relation; 

NSEG  *  maximum  number  of  history  segments  associated  with  a  single 
prescribed  loading  history; 

NINC  »  maximum  number  of  increments  into  which  a  single  prescribed 
loading  history  has  been  subdivided; 

NVAR  >  number  of  different  response  parameters  (such  as  q,  p',  u, 
Ej  and  ey)  used  to  define  the  results  of  the  various 
experimental  tests;  and, 

NRLN  »  number  of  different  response  relations  (such  as  q  vs  p",  q 
vs  p'  vs  €j,  u  vs  Gj  and  ey  vs  e^)  used  to  define  the 
results  of  the  various  experimental  tests. 

The  program  is  currently  set  up  to  accommodate  five  different  response 
relations  (NRLN  =  5),  which,  as  a  group,  are  fully  defined  by  means  of  five 
distinct  response  parameters  (NVAR  =  5).  If,  for  some  reason,  it  becomes 
necessary  to  use  relations  and/or  parameters  that  are  different  from  those  listed 
above  to  characterize  the  results  of  a  given  experimental  test,  minor  changes 
must  be  made  to  the  coding  which  handles  the  storage  and  retrieval  of  the 
model  predictions.  These  coding  changes  are  quite  straightforward,  and  are 
restricted  only  to  subroutines  EVAL,  PLTCHK  and  RSDUAL.  The  remaining 
dimensions  are  currently  set  at  NDIM  =  17,  NTST  =  6,  NPLT  =  3,  NEXP  =  20, 
NSEG  =  4  and  NINC  =  60. 
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XVQ.  EXAMPLE  1:  MODEL  CALIBRATION 

To  provide  a  first  example  of  the  code's  capabilities,  MODCAL  has  been 
used  to  assist  in  the  calibration  of  the  Bounding  Surface  plasticity  model  for  a 
particular  cohesive  soil.  As  may  be  seen  from  the  input  and  output  files  which 
follow,  the  program  was  run  to  establish  the  optimal  values  of  the  five  model 
parameters  G,  Rc,  Ac,  c  and  hc  for  the  case  of  the  laboratory  prepared  Kaolin 
employed  by  Jafroudi  (1983)  (see  also  Herrmann  et  al  (1981b)).  The  specified 
calibration  data  base  includes  the  results  of  conventional  undrained  triaxial 
compression  tests  (in  the  form  of  the  experimentally  observed  q  vs  e^,  p'  vs 
and  u  vs  responses)  performed  on  specimens  at  overconsolidation  ratios  of 
OCR  =  1,  2  and  6.  All  experimental  observations  are  weighted  equally,  and  the 
objective  function  is  formed  by  using  the  absolute-Euclidean  measure  of  error 
(KNRM=0  and  KRSL=0),  together  with  the  values  PLIM  =  0.01  and  TLIM  =  0.20. 


S . 

EXAMPLE  i:  MODEL  CALIBRATION.  JAFROUPI  (1983)  KAOLIN 

--  0CR=li2.6. 

v  3  1 

M  0.130 

0  0 
0.018 

0  0 
1.180 

5  5 

0.737 

9  10  14  16 

5000.00  1000000.0 

101.35 

101.35 

2.400 

0.100 

0.050 

1.000 

1.000 

0.400 

1.000 

0.500 

i 

•  n 

i  l  .0 

0.500 

0.68 

0.69 

0.72 

0 

1.000 

392.20 

392.20 

392.20 

1.00 

1.000 

392.20 

196.10 

65.37 

0.01 

1.00 

2.00 

6.00 

0.01 

0.20 

r  i  2 

0  392.2  0 

392.2  1 

0.02  1 

1 

1 

20 

1.00 

Lv  0  392.2  0 

392.2  1 

0.12  1 

1 

1 

20 

1.00 

2  2 

0  196.1  0 

196.1  1 

0.02  1 

1 

1 

20 

1.00 

V  0  196.1  0 

196.1  1 

0.12  1 

1 

1 

20 

1.00 

3 

0  65.4  0 

65.4  1 

0.02  1 

1 

1 

20 

1.00 

0  65.4  0 

65.4  1 

0.12  1 

1 

1 

20 

1.00 

1  3 

2  15 

0.0 

65.0 

3  15 

105.0 

155.0 

4  15 

176.0 

187.0 

202.0 

208.0 

f  214.0 

216.0 

217.0 

217.0 

215.0 

212.0 

206.0 

0.00 

0.12 

0.25 

0.50 

0.75 

1.00 

1.50 

2.00 

3.00 

4.00 

5.00 

6.00 

8.00 

10.00 

12.00 

392.2 

382.9 

362.2 

329.9 

304.9 

284.5 

256.5 

236.5 

216.5 

204.2 

196.5 

191.5 

194.9 

179.9 

174.9 

■  0.00 

^  3.00 

0.12 

0.25 

0.50 

0.75 

1.00 

1.50 

2.00 

4.00 

5.00 

6.00 

8.00 

10.00 

12.00 

0.0 

30.0 

65.0 

114.0 

146.0 

170.0 

203.0 

225.) 

Si  247.0 

260.0 

268.0 

273.0 

279.0 

283.0 

286.0 

S  0.00 

0.12 

0.25 

0.50 

0.75 

1.00 

1.50 

2.00 

3.00 

i  2  3 

*  2  15 

0.0 

4.00 

5.00 

6.00 

8.00 

10.00 

12.00 

50.0 

3  15 

70.0 

135.0 

4  15 

155.0 

169.0 

180.0 

187.0 

v  195.0 

199.0 

202.0 

202.0 

202.0 

200.0 

198.0 

0.00 

0.12 

0.25 

0.50 

0.75 
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EXAMPLE  i:  MODEL  CALIBRATION.  JAFROUDI  (1983)  KAOLIN  --  0CR=1>2>6. 


INPUT  data: 
CONTROL  CODES*. 


JOPT  =  1 
JRUN  =  0 
JPLT  =  0 
NTST  =  3 
KRSL  =  0 
KNRM  =  0 
POUR  =  l.OOE+OO 
PLIM  =  1.09E-02 
TLIM  =  2.00E-01 
NOPT  *  5 
KOPT  5  5  9 


input  data: 

MATERIAL  PROPERTIES*. 


LAMBDA 
KAPPA 

NUC  <CSL  SLOPE) 

MUE/MUC 

SHEAR  MODULUS  G 
CBULK  MODULUS  GAMMA  «  1.000E+06 
TRANSITIONAL  STPESS  PL  *  1.013E+02 
ATMOSPHERIC  PRESSURE  PA  -  1.013E+02 


MODEL  constants: 


RC  =  2.400E+00 
AC  *  l.OOOE-Ol 
T  *  5.000E-92 
RE/RC  *  1.000E*00 
AE/AC  *  1.000E400 
PROJECTION  POINT  C  *  4.000E-01 
ELASTIC  NUCLEUS  S  =  l.OOOE+OO 
HC  =  5.000E-01 
MC  =  5.000E-01 
HE/HC  *  l.OOOE+OO 
ME/MC  *  l.OOOE+OO 


10  14  10 


*  1.300E-01 

*  1.800F-02 

*  1.180E+00 
«  7.370E-01 
=  5.000E+03 


INITIAL  STATE  PROPERTIES 


i.'. 


I 


U 


ti 


♦ 


li 


TEST 

OCR 

V01D-R 

P-PRECON 

P-CONFIN 

1 

1.00 

0.68 

3.922E+02 

3.922E+02 

2 

2.00 

0.69 

3.922EF02 

1.961E402 

3 

6.00 

0.72 

3.922E402 

6.537E401 

INPUT  data: 

ITERATION  INFORMATION  AND  ANALYTICAL  OPTIONS*. 


HAX  #  OF  ITERATIONS  -  10 

MAX  RELATIVE  ERRORS  -  1.000E-02 

ACCELERATION  FACTOR  LIMITS  *  l.OOOE+OO 

OtS  UNTRAINED  CONDITIONS  Sttt 

SIS*  ENGINEERING  STRESSES  AND  STRAINS  ASSUKD  Sttt 

mt  REFORMULATED  NEARLY  INCOMPRESSIBLE  ANALYSIS  SSSS 


INPUT  data: 

SPECIFIED  LOADING  HISTORIES*. 


TEST 

SEG-# 

♦-INC 

S-RATIO 

SS-CODE 

SIG/EPS-X 

SIG/EPS-Y 

SI6/EPS-Z 

TAU/GAM-XY 

TAU/GAH-XY 

1 

1 

20 

1.00 

001111 

3.922E402 

3.922E402 

2.000E-02 

O.OOOE-Ol 

O.OOOE-OI 

1 

2 

20 

1.00 

001111 

3.922E402 

3.922E402 

1.200E-01 

O.OOOE-Ol 

O.OOOE-Ol 

2 

1 

20 

1.00 

001111 

1.961E402 

1.961E402 

2.000E-02 

O.OOOE-Ol 

O.OOOE-Ol 

2 

2 

20 

1.00 

001111 

1.961E402 

1.961E402 

1.200r-01 

O.OOOE-Ol 

O.OOOE-Ol 

3 

1 

20 

1.00 

001111 

6.540E401 

6.540E401 

2.000E-02 

O.OOOE-Ol 

O.OOOE-Ol 

3 

2 

20 

1.00 

001111 

6.540EF01 

6.540E401 

1.200E-01 

O.OOOE-Ol 

O.OOOE-Ol 
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TAU/GAH-YZ 


O.OOOE-Ol 

O.OOOE-Ol 

O.OOOE-Ol 

O.OOOE-Ol 

O.OOOE-Ol 

O.OOOE-Ol 


INPUT  data: 

SPECIFIED  EXPERIMENTAL  OBSERVATIONS  AND  WEIGHT INS  FACTORS 


INPUT  mta: 

CALIBRATION  CONTROL  CODES: 


NUMBER  OF  DIMENSIONS  -  3 
MAXIMUM  NUMBER  OF  FUNCTION  EVALUATIONS  *  400 
OUTPUT  DATA  PRINT  CODE  =  2 
C0NVER6ENCE  CRITERION  CODE  =  3 
CUBIC  INTERPOLATION  OPTION  CODE  =  0 
HESSIAN  CONDITIONING  OPTION  CODE  •  0 
•  OF  NON-ZERO  HESSIAN  VALUES  TO  BE  READ  =  0 
PARAMETER  SCALING  OPTION  CODE  =  0 


«  OF  IMPROVED  STARTING  ESTIMATES  REQUIRED  =  1 

ABSOLUTE  CONVERGENCE  CRITERION  FOR  F(X)  *  1.000E-02 

ABSOLUTE  CONVERGENCE  CRITERION  FOR  SIX)  *  1.000E-08 

VALUE  OF  THE  LINE-SEARCH  PARAMETER!  SIGMA  •  9.000E-01 
VALUE  OF  THE  BRACKET-CHECK  PARAMETER*  ROA  *  1.000E-04 
VALUE  OF  THE  BRACKET-CHECK  PARAMETER.  TAU  -  1.000E-01 
VALUE  OF  THE  LINE  SEARCH  EXIT  CRITERION  =  1.000E-03 

VALUE  OF  THE  FAIL-SAFE  GLOBAL  CRITERION  *  1.000E-04 

I  OF  FN-EVALUATIONS  BEFORE  FAIL-SAFE  EXIT  *  100 
ERROR  BOUND!  MACHINE  REPRESENTATION  OF  XX  =  1.000E-09 
ERROR  BOUND:  MACHINE  REPRESENTATION  OF  FX  -  1.000E-04 
MINIMUM  VALUE  OF  THE  DIFFERENCE  INTERVALS  =  5.000E-03 
MAXIMUM  ERROR  IN  FUG-DIFFERENCE  CALCULATN  *  1.000E-04 
SCALING  FACTOR  FOR  THE  OBJECTIVE  FUNCTION  •  l.OOOE+OO 


INPUT  data: 

INITIAL  VALUES  OF  THE  FUNCTION  VARIABLES  . . . 


5.00E+03  2.55E+00  1.00E-01  4.00E-0J  5.00E-01 
INPUT  data: 

LOUER  AND  UPPER  BOUNDS  ON  THE  FUNCTION  VARIABLES  X1.X2.....XN5 


2.50E+03  2.00E400  3.00E-02  0.00E-01  5.00E-02 

1.00E+04  3.00E400  2.00E-01  8.00E-01  2.00E+00 

INPUT  data: 

INITIAL  VALUES  OF  THE  DIFFERENCING  INTERVALS  D1»D2*,...DN: 
5.00E402  5.00E-02  1.00E-02  3.00E-02  5.00E-02 
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OUTPUT  mta: 

THE  RESULTS  Of  THE  OPTIMIZATION  ROUTINE  FOLLOW: 


l-I  l-F  l-G 

FI 

F2  SI  S2 

XI 

X2 

X3 

X4 

0  6 

1 

0.000E-01 

1.209E-01  0.000E-01  -1.162E-01 

5.000E4Q3 

2.551E+00 

1.000E-01 

4.000E-01 

1  17 

2 

1.209E-01 

1.140E-01  -1.162E-01  3.986E-02 

5.918E+03 

2.520E400 

7.425E-02 

4.984E-01 

1.790E-0: 

2  28 

3 

1.140E-01 

1.057E-01  -2.221E-02  -3.672E-03 

6.006E+03 

2.534E+00 

6.624E-02 

4.407E-01 

2.633E-0: 

3  44 

4 

1.057E-01 

1.057E-01  -8.100E-03  -8.479E-03 

6.006E+03 

2.534E400 

6.624E-02 

4.407E-01 

2.633E-0 

4  64 

5 

1.057E-01 

1.057E-01  -7.830E-01  9.880E-01 

6.006E+03 

2.534E400 

6.624E-02 

4.407E-01 

2.633E-0 

6  112 

8 

1.020E-01 

I.020E-0!  -1.938E400  -6.981E-01 

6.475E+03 

2.574E+00 

6.393E-02 

4.441E-01 

3.008E-0 

7  126 

9 

1.020E-01 

1.020E-01  -5.955E-04  3.318E-03 

6.476E+03 

2.574E+00 

6.393E-02 

4.441E-01 

3.008E-0 

8  137 

10 

1.020E-01 

1.020E-01  -1.210E-06  2.506E-05 

6.474E+03 

2.574E400 

6.392E-02 

4.441E-01 

3.009E-0 

9  148 

11 

1.020E-01 

1.020E-01  -2.323E-05  7.391E-06 

6.475E403 

2.574E+00 

6.393E-02 

4.441E-01 

3.008E-0 

Util  THE  INVERSE  HESSIAN  HAS  BEEN  RESET  TO  THE  IDENTITY  MATRIX  ttttt 

10  171 

12 

1.020E-01 

1.008E-01  -9.659E-02  2.457E-02 

6.519E+03 

2.571E+00 

6.243E-02 

4.587E-01 

2.607E-0 

11  184 

13 

1.008E-01 

1.008E-01  -3.280E-02  -5.474E-03 

6.546E+03 

2.570E40C 

5.712E-02 

4.656E-01 

2.554E -0 

12  199 

14 

1.008E-01 

1.008E-01  -8.714E-04  -3.526E-04 

6.545E+03 

2.570E+00 

5.711E-02 

4.655E-01 

2.555E-0 

13  214 

15 

1.008E-01 

1.008E-01  -9.221E-04  -1.240E-03 

6.545E403 

2.570E+00 

5.711E-02 

4.656E-01 

2.556E-0 

14  227 

16 

1.008E-01 

1.008E-01  -2.456E-04  -1.185E-04 

6.545E+03 

2.570E+00 

5.711E-02 

4.656E-01 

2.554E-0 

15  241 

17 

1.008E-01 

1.008E-01  -2.144E-04  -3.1B8E-04 

6.545E+03 

2.570E+00 

5.712E-02 

4.656E-01 

2.556E-0 

16  267 

19 

1.007E-01 

1.009E-01  -5.078E-04  -5.07BE-04 

6.542E+03 

2.573E+00 

5.732E-02 

4.651E-01 

2.576E-0 

ttttt  THE  SEARCH  HAS  BEEN  TERMINATED  DUE  TO  INSUFFICIENT  PROGRESS  ttttt 


xvra.  EXAMPLE  m  MODEL  PREDICTION 

To  provide  a  second  example  of  the  code's  capabilities,  MODCAL  has 
been  used  to  generate  model  predictions  for  a  typical  set  of  experimental  test 
conditions.  As  may  be  seen  from  the  input  and  output  files  which  follow, 
predictions  are  required  for  three  distinct  loading  histories.  These  three  loading 
histories  simulate  the  conventional  undrained  triaxial  compression  tests  performed 
by  Jafroudi  (1983)  on  samples  of  Kaolin  at  overconsolidation  ratios  of  OCR  = 
1,  2  and  6  (see  also  Herrmann  et  al  (1981b)).  For  completeness,  the  relation 
plotting  option  has  also  been  invoked  (3PLT  =  1),  and  in  the  nine  plots  which 
result,  the  symbol  represents  an  experimental  observation,  while  the  symbol 
"•*"  denotes  a  model  prediction. 


EXAMPLE  II:  MODEL  PREDICTION.  JAFROUD!  (1983)  KAOLIN  --  OCR^l»2.6. 
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INITIAL  STATE  PPCPFPTIES: 
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UNCLASSIFIED 


USER'S  MANUAL  FOR  MODCAL  -  BOUNDING  SURFACE  SOIL 
PLASTICITV  MODEL  CALI BRA.  .  (U>  CALIFORNIA  UNIV  DAVIS 
DEPT  OF  CIVIL  ENGINEERING  J  S  DENATALE  ET  AL.  FEB  83 
NCEL-CR-82.  811  N62474-82-C-8276  .  F/G  8/13 


RESOLUTION  TEST  CHART 

JREAU  OF  STANDARDS-  1963-A 
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o  n  n  nnnnnnnnn  non  nnnnnnnnnnnnnnn 


MODCAL 


THIS  SUBROUTINE  SERVES  AS  THE  MAIH  DRIVING 
PROGRAM  FOR  THE  BOUNDING  SURFACE  SOIL  PLASTICITY 
MODEL  CALIBRATION  AND  PREDICTION  CODE. 

WRITTEN  BY  J.S.  DE  NATALE, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

VERSION  I:  MAY  1982. 

DIMENSION  XV (17) 

nt  t  npeu 

CALL  DATAIN ( JOPT , JRUN ) 

IF(JOPT.NE.O)  CALL  BSMOPT(XV) 

IF(JRUN.NE.O)  CALL  EVAL(XV,FV, 1 ) 

CALL  EXIT 
END 


SUBROUTINE  DATAIN(IOPT,IRUN) 

THIS  SUBROUTINE  READS  IN  ALL  VALUES  REQUIRED  TO  CALIBRATE 
AND/OR  USE  THE  BOUNDING  SURFACE  PLASTICITY  MODEL  PROGRAMS. 

WRITTEN  BY  J.S.  DE  NATALE, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

VERSION  Is  MAY  1982. 

COWWN/BLX7/SFUN ,  SF  ( 1 7 ) ,  XL  ( 1 7 ) ,  XU  ( 1 7 ) 

COMMON /BLXA/NTST , JOPT , JPLT , NOPT , XOPT (17) 

COMMON/BLKB/W1 (6),W2(6, 3)>V3(6, 3»20) 

COMMON/BLKC/PRP 1 ( 1 9 ) , PRP2 ( 6 , 4 ) , IDAT ( 3 ) » RDAT (2 ) , 

•  NSEG ( 6 ) , LTYP (6,4,7), VALU (6,4,7) 
COMMON/BLKD/NPLT (6 ) , KPLT ( 6 , 3 ) , NEXP ( 6 , 3 ) » 

•  XV(6,3,20),YV(6,3,20),PINC(61,5) 

COMMON /BLKE/KRSL , POWR , XMXR (5 ) , YMXR(5), XMXP (6, 3 ) , 

•  YMXP(6, 3)»PLIM,TLIM 
DIMENSION  TITL(20) 

FORMAT  STATEMENTS 

800  FORMAT (20A4) 

804  FORMAT (1615) 

808  FORMAT (8E 10. 3) 

812  FORMAT(1I5,5X,4E10. 3) 

816  FORMAT(3I5,5X,4E10.3) 

820  FORMAT (4 (2I5,E10.3)) 
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f 


1 


824  FORMAT(6(I1,E9.2),5X,I5,E10.3) 
828  F0RMAT(4 (213,14,810. 3)) 

900  FORMAT (1H1/5X.20A4) 

904  FORMAT(//5X, 

5X, 

15X, 

15X, 

15X, 

15X, 

15X, 

15X, 

15X, 

15X, 

15X, 

15X, 

15X, 

908  FORMAT (//5X, 

5X, 


912 


916 

920 


15X, 
15X, 
15X, 
15X, 
15X, 
15X, 
FORMAT (//5X, 
5X, 
15X, 
*X, 
2(4X 


FORMAT (//5X, 
’AND 
5X, 
15X, 
15X, 
15X, 

924  F0RMAT(/10X, 
928  FORMAT (/I OX, 
932  FORMAT ( / 1 0X , 
936  FORMAT(/10X, 
940  F0RMAT(/10X, 


=  ’,13/ 
=  '  ,13/ 
=’,13/ 
=  ',13/ 


INPUT  DATA: */5X, ’CONTROL  CODES:'/ 

- ’// 

JOPT 
JRON 
JPLT 
NTST 
KRSL 
KNRM 
POUR 
PLIM 
TLIM 
NOPT 
KOPT 


s’, 13/ 

13/ 

s', 1PE10.2/ 
s’,  810.2/ 

s’,  810.2/ 

=',13/ 

=’, 15(13, 3X>) 

INPUT  DATA: */5X, ’MATERIAL  PROPERTIES:'/ 
- •  // 


15X, 'LAMBDA 

1PE10.3/ 

15X, ’KAPPA 

s’, 

E10.3/ 

15X, 'MUC  (CSL  SLOPE) 

s’, 

E10.3/ 

15X, 'HUE /MUC 

s’. 

E10.3/ 

15X, 'SHEAR  MODULUS  G 

-  , 

810.3/ 

15X,'CBULK  MODULUS  GAMMA 

s’, 

E10.3/ 

15X, ’TRANSITIONAL  STRESS  PL 

- 1 

BIO. 3/ 

15X, 'ATMOSPHERIC  PRESSURE  PA 

s’, 

E10.3/// 

51 , ’MODEL  CONSTANTS : ’/5X, 
15X,'RC  =' ,E10. 3/15X, 'AC 

mmmmm 

s’,E10.3/ 

»// 


T  s' , BIO. 3/15X, 'RE/RC  s»,E10-3/ 
AE/AC  s’ ,E10. 3/ 

PROJECTION  POINT  C  s’, El 0.3/ 

ELASTIC  NUCLEUS  S  s’, El 0.3/ 

HC  s' ,E10. 3/15X, ’MC  s*,E10.3/ 
HE/HC  s' ,E10. 3/15X, 'ME/MC  s’,E10.3) 
INITIAL  STATE  PROPERTIES:’/ 

- ’// 

TEST' ,5X, ’OCR’ , 2X, ’VOID-R’ ,4X, ’P-PRECON’ , 

P-CONFIN’ /15X, ' - ’ , 5X, ' - ',2X, ' - ', 

’ - ’)/) 


FORMAT (15X,l4,2F8. 2, 1P3E12.3) 


INPUT  DATA :’/5X, 'ITERATION  INFORMATION 
ANALYTICAL  OPTIONS:'/ 


'// 


MAX  #  OF  ITERATIONS 
MAX  RELATIVE  ERRORS 
ACCELERATION  FACTOR 


s’, 13/ 
s’ , 1PE10. 3/ 
LIMITS  s’,B10.3) 


••••  DRAINED  CONDITIONS  •••••) 

••••  UNDRAINED  CONDITIONS  •••••) 

•••*  engineering  stresses  and  STRAINS  ASSUMED  ••••’) 
••••  true  STRESSES  AND  NATURAL  STRAINS  ASSUMED  ••••’) 
••••  REFORMULATED  NEARLY  INCOMPRESSIBLE  *, 

•  'ANALYSIS  ••••’) 

944  FORMAT(/10X, '••*•  NON-REFORMULATED  NEARLY  INCOMPRESSIBLE  ’, 

•  ’ANALYSIS  ••••») 
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948  FORMAT (//5X, ’INPUT  DATA: */5X, ’SPECIFIED  LOADING  HISTORIES:*/ 


•  51,  • - •// 

•  5X, *TEST\2X, »SEG-#' ,2X, *#-INC' ,2X, 'S-RATIO* , 

•  2X, 'SS-CODE' , 3X, *SIG/EPS-X* , 3X, ’SIG/EPS-T* , 

•  3X, 'SIG/EPS-Z  • ,  2X, *TA0/GAM-XY*,2X, 'TAU/GAM-XY* , 

•  2X, ' TAU/GAM-YZ 1  /5X , ' - *,2(2X, ' - •), 

•  2(2X  •———•)  3(3x  3(2X 

952  FORMAT (5X , I*, 2 (2X| 15 ),F9! 2,3x1611, 1P6E12.3) 

956  FORMAT (20X) 

READ  HEADING  INFORMATION  AND  CONTROL  CODES 
READ (5, 800)  TITL 

RE AD ( 5 , 804 )  NTST , JOPT , JRUN , JPLT , KRSL , KNRM, NOPT , 

•  (KOPT (II), 11=1, NOPT) 

IOPT=JOPT 

IRONs JRUN 

IF (KNRM . EQ . 0 )  POUR=1.00 
IF(KNRM.NE.O)  POVR=2.00 

READ  MATERIAL  AND  MODEL  PROPERTIES 

READ (5, 808)  (PRP1 (II), 11=  1,15) 

READ (5, 808)  (PRP1 (II),II=16, 19) 

IF(PRP1(1*)-GT. 1.00)  PRP1 ( 14 )=1. 00 
IF(PRP1 (15).LT. 1. 00)  PRP1 (15)*1.00 
DO  110  IIsl.NTST 

110  READ (5, 812)  ITST, (PRP2(ITST, JJ),JJs1,4) 

READ  CONVERGENCE  CRITERIA  AND  ITERATION  INFORMATION 

READ(5,816)  (IDAT (II ) , 11= 1 , 3 ) , (RDAT (II ) , 11= 1,2), PLIM,TLIM 

READ  SPECIFIED  LOADING  HISTORIES 

DO  120  11=1, NTST 

READ (5, 804)  ITST,NSEG(II) 

DO  120  JJ=1,NSFG(II) 

READ (5, 824 )  (LTYP(II,JJ,nO,VALU(II,JJ,!CK),KKr1,7) 

120  IF(VALO(II,JJ,7).EQ.O.O)  VALO(II, JJ,7)=1. 0 

READ  EXPERIMENTAL  DATA  AND  PLOTTING  INSTRUCTIONS 

HSNTsO 

DO  140  11=1, NTST 

READ (5, 820)  ITST,NPLT(II ),W1 (II),NWPT 

READ (5, 820)  ( (KPLT(II,KK),NEXP(II,KK),W2(II,KK)), 

•  KK=1,NPLT(II)) 

NSWT=NSWT+NNPT 

DO  130  JJ=1,NPLT(II) 

IF(NEXP(II , JJ ) . EQ. 0 )  QO  TO  140 

READ (5, 808)  (YV(II, JJ,RR),KR=1,NEXP(II, JJ)) 

READ (5, 808)  (XV(II,JJ,KK),KK=1,NEXP(II,JJ)) 

130  CONTINUE 


140  CONTINUE 

DO  150  11=1, NTST 

IF(H1 (II).BQ.O. 0)  N1(II)=1.0 

DO  150  JJ=1,NPLT(II) 

IF(H2(II,JJ).EQ.0.0)  W2(IIfJJ)=1.0 
DO  150  KK=1,NEXP(II,JJ) 

150  H3(II,JJ,KX)s1-0 
C 

C  READ  SPECIAL  EXPERIMENTAL  WEIGHTINGS 
C 

IF(NSWT.EQ.O)  GO  TO  170 

A1=FLOAT(NSVT)/4.0 

A2=AINT(A1) 

IFU1-A2.GT.0.0)  A2=A2+1.0 
NW=IFIX(A2) 

NX=NSWT-4«(NW-1) 

DO  160  11=1, NV 

IF(II.LT.NW)  READ (5 ,828)  ((IT, IP, IE, 

•  W3(IT,IP,IB)),JJ=1,  4) 

160  IF(II.BQ.NW)  READ (5, 828}  ((IT, IP, IE, 

•  W3 (IT, IP, IB) ),JJ=1 ,NX) 

170  CONTINUE 

C 

C  INITIALIZE  SCALING  FACTORS 
C 

DO  180  11=1, NOPT 
180  SF(II)=1.0 
C 

C  ASSIGN  DEFAULT  VALUES 

C 

IF(PLIM.LT.O.OI)  PLIMsO.OI 
IF(TLIM.LT.O.OI)  TLIM=0.01 
IF(PLIM.GT. 1.00)  PLIM= 1.00 
IFOTLIM.GT. 1.00)  TLIM=1.00 
IF(IDAT(2).LE.  0)  IDAT(2)=10 
IF(IDAT(2).GT.  20)  IDAT(2)=20 
IF(RDAT(1).LE.0.0)  RDAT(1)=0.30 
IF (RDAT ( 2 ) . LE • 0. 0 )  RDAT(2)=0.01 
C 

C  PRINT  INPUT  DATA 

C 

WRITE (6, 900)  TITL 

WRITE ( 6 , 904 )  JOPT , JRUN , JPLT , NTST , KRSL , RNRM , 

•  POWR , PLIM , TLIM , NOPT , (XOPT (II), 11*1, NOPT) 
WRITE (6, 908)  (PRP1 (II), 11=1, 19) 

WRITE(6,912) 

DO  190  11=1, NTST 

190  WRITE(6,916)  II,PRP2(II,4 ), (PRP2(II, JJ ) , JJ= 1 , 3) 

WRITE (6, 920)  IDAT(2),RDAT(2),RDAT(1 ) 
IF(PRP1(6).EQ.O.O)  WRITE(6,924) 

IF(PRP1 (6).NE. 0. 0)  WRITE (6, 928) 

IF(IDAT(3).EQ.  0)  WRITE(6,932) 

IF(IDAT(3).BQ.  1)  WRITE(6,936) 

IF(IDAT(1 ).EQ.  0)  WRITE (6, 940) 


o  o  o  o  n  o  oonoo  ooo  o  o  o  o 


IF(IDAT(1 ).EQ.  1)  WRITE (6, 944) 

WRITE (6, 948) 

DO  210  11*1, NTST 
DO  200  JJ*1,NSEG(II) 

WRITE (6, 952)  II , J J , LTYP (II , J J , 7 ) , VALU (II ,  J  J ,  7 ) , 

•  (LTYP (II , J J , XX ) , XX* 1,6), 

•  (7ALU(II,JJ,XX),EX=1,6) 

200  CONTINUE 

WRITE (6 ,956) 

210  CONTINUE 

ESTABLISH  EXPERIMENTAL  MAXIMA 
IF  MODEL  CALIBRATION  IS  REQUIRED 

IF(JOPT.EQ.I)  CALL  FNDMAX (NTST , NSWT ) 

RETURN 

END 


SUBROUTINE  FNDMAX (NT, NS) 

THIS  SUBROUTINE  ESTABLISHES  THE  MAXIMA 
OF  THE  SPECIFIED  EXPERIMENTAL  OBSERVATIONS, 

AS  REQUIRED  BY  SUBROUTINE  DAT AIN. 

COMMON /BLKB/W 1 ( 6 ) , W2 ( 6 , 3 ) # W3 ( 6 , 3 , 20 ) 

COMMON /BLKD/NPLT (6) , KPLT(6 , 3 ) , NBXP(6, 3 ) , 

•  XV(6,3,20),TV(6,3,20),PINC(61,5) 
COMMON/BLKE/KRSL , POWR , XMXR ( 5 ) , YMXR (5 ) , XMXF ( 6, 3 ) , 

•  YMXP(6, 3),PLIM,TLIM 

ESTABLISH  MAXIMA 

DO  100  11*1,5 
XMXR (II) =0.0 
100  YMXR(II )*0. 0 
DO  120  11*1, NT 
DO  120  JJ=1,NPLT(II) 

KR*KPLT(II, JJ ) 

XMXP(II , JJ )*ABS(XV (II, JJ , 1 ) ) 

YMXP(II , JJ )*ABS (YV (II , JJ , 1 ) ) 

DO  110  KK*2,NEXP(II,JJ) 

IF(XMXP(II , JJ ). LT. ABS (XV (II , JJ ,KK) ) ) 

•  XMXP(II.;j)*ABS(XVdI,JJ,EK)) 

110  IF(YMXP(IX, JJ  LT.ABSC"  II,JJ,XX))) 

•  YMXP(II,..y  ABf  V(II,JJ,XX)) 

IF(XMXR  (KR  ).LT.  jlt_aP(II,  JJ  ) )  XMXR  (KR  )=XMXP  (II ,  J  J  ) 

120  IF(YMXR(RR).LT.YMXP (II, JJ))  YMXR(KR)*YMXP(II, JJ) 

PRINT  EXPERIMENTAL  MAXIMA 

WRITE (6, 1 30 ) 

130  FORMAT (//5X, 'INPUT  DATA: '/5X, 'SPECIFIED  EXPERIMENTAL  ', 
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o  o  o  o  o  o 


•  'OBSERVATIONS  AND  WEIGHTING  FACTORS:'/ 

•  5X, * - 

•  ' - «// 

•  5X, ' TEST ' , 31 , * PLOT » , 3X , ' KIND ' , 2X , ' WF-TEST ’ , 

•  2X , ' WF-PLOT ' , 2X, ' #-PTS ' , 2X, ' IMAX (PLOT ) • , 

•  21, 'XMAX (PLOT ) ' , 2X, 'IMAX (KIND) ', 2X, 'XHAX(KIND)'/ 

•  51, ' - '  ,2(3X, ' - •), 

•  2(2X,  ’ - ’  ),2X, ' - \4(2X, ' - •)/) 

DO  160  11*1 , NT 

DO  140  JJ*1 ,NPLT(II ) 

KR=KPLT(II,JJ) 

140  WRITE (6, 150)  II, JJ,KR,W1 (II),W2(II, JJ),NEXP(II, JJ), 

•  YMXP (II , J J ) , XMXP (II , J J ) , TMXR (KR ) , XMXR (KR ) 
150  F0RHAT(2X,3(3X,I4),2F9.3.2X,I5,1P4B12.3) 

160  WRITE(6,170) 

170  FORMAT (20X) 

PRINT  SPECIAL  EXPERIMENTAL  WEIGHTINGS 

IF(NS.EQ.O)  GO  TO  210 
IS=0 

WRITE(6, 180) 

180  FORMAT (//5X, 'INPUT  DATA: '/5X, 'SPECIAL  EXPERIMENTAL  ', 

•  'WEIGHTINGS: '/5X, 32 ('-’ )//5X, 'PT-#' ,3X, 'TEST' , 

•  3X, 'PLOT* ,2X, 'POINT' ,6X, 'WEIGHT' /2X, 3 (3X, • - ' ), 

•  21,  • - ' , 6X , ' - •/) 

DO  190  11*1 ,NT 

DO  190  JJsl.NPLT(II) 

DO  190  KK*1,NEXP(II,JJ) 

WT=W3(II,JJ,KK) 

IF(WT.NE.I.O)  IS*IS+1 

190  IF(WT.NE.I.O)  WRITE (6, 200)  IS, II, JJ.KK.WT 
200  F0RMAT(2X,4 (4X,I3), 1 PE 12. 3) 

210  CONTINUE 
RETURN 
END 


SUBROUTINE  EVAL(XOPT,RSDL, JCOD) 

C 

C  THIS  SUBROUTINE  PERFORMS  A  SINGLE-ELEMENT, 

C  INCREMENTAL-ITERATIVE  FINITE  ELEMENT  ANALYSIS 

C  FOR  HOMOGENEOUS  LOADING  CONDITIONS.  FOR  USE 

C  WITH  MATERIAL  MODELS  SUCH  AS  THE  BOUNDING 
C  SURFACE  SOIL  PLASTICITY  MODEL. 

C 

C  WRITTEN  BY  L.R.  HERRMANN, 

C  RECODED  BY  J.S.  DE  NATALE, 

C  DEPARTMENT  OF  CIVIL  ENGINEERING, 

C  UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

C  LAST  REVISED:  MAY  1982. 

C 

COMMON /BLK7 /SFUN , SF ( 1 7 ) , XL ( 1 7 ) , XU ( 1 7 ) 


o  n  o  o  o  o  o  o  o 


COMMON /BLK A /NTST , JOPT , JPLT , NOPT , KOPT (17) 
COMMON/BLKB/V1(6),W2(6,3),W3(6,3,20) 

COMMON/BLKC/PRP 1 ( 19 ) , PRP2 ( 6 , 4 ) , IDAT ( 3 ) , RDAT (2 ) , 

•  NSEG ( 6 ) , LTYP ( 6 , 4 , 7 ) » VALU (6,4,7) 

COMMON /BLKD /NPLT ( 6 ) , KPLT ( 6 , 3 ) , NEXP ( 6 , 3 ) , 

•  XV(6,3,20),YV(6,3,20),PINC(61,5) 

COMMON /BLKE/KRSL , POVR , XMXR ( 5 ) , YMZR ( 5 ) , XMXP (6,3), 

•  YMXP(6,3),PLIM,TLIM 

DIMENSION  ILOD ( 6 ) , VLOD ( 6 ) , DLOD ( 6 ) , SIGB ( 6 ) , EPSB ( 6 ) , DSIG ( 6 ) , 

•  DEPS(6) ,STOR(6),PROP(21 ),C(6,6),S(7,7)»R(7),RP(6) 
DIMENSION  XOPT (17) 

FORMAT  STATEMENTS 

900  FORMAT ( 1H 1, 2X, 'OUTPUT  DATA: '/3X, 'MODEL  PREDICTIONS  ', 

•  'CORRESPONDING  TO  TEST  DATA  SET  #»,I2, '  ARE  AS  POLLOWS: 

•  /3X.17C - ') 

•  ///3X, 'N' ,4X, 'EPS-X* ,4X, 'EPS-Y' ,4X, 'EPS-Z' , 3X, 'GAM-XY' , 

•  3X , ' GAM-XZ ' , 3X , ' GAM- YZ • , 5X , ' SIG-X ' , 5X , ' SIG-Y ' , 5X , ' SIG-Z 

•  4X , • TAO-XY • , 4X , ' TAU-XZ • , 4X , ' TAU-YZ • , 6X , *  0 • , 3X , ' #-IT • / 

•  3X, ,  3(41, ' - '  ),3(3X, ' - ' ),  3(51, ' - '), 

•  3(4X, ' - '  ),5X, ' - '  ,21, ' - •/) 

904  FORMAT* //3X, ’CONVERGENCE  (ID  NOT  OCCUR  FOR  INCREMENT* ,13, 

•  •  -  ERSIG=\E9.3,'  AND  EREPS=' ,E9. 3//) 

908  F0RMATOX.I3, 1P6E9. 1.7E10.2,  IX, 13) 

912  FORMAT ( / ) 

916  FORMAT (1H1,2X, 'OUTPUT  DATA: ' /3X, ’MODEL  PREDICTIONS  ’, 

•  '—IN  TERMS  OF  THE  KEY  RESPONSE */3X, 'PARAMETERS  —  ', 

•  'CORRESPONDING  TO  TEST  DATA  SET  #',I2, '  ARE  AS  FOLLOWS: 

•  /3X.16C - ')//l6X, 'N',5X, 'EPS-1 'f9X,’Q*,9X,'P',9X,'U' 

•  5X, 'E-VOL' /16X, '-' ,5(5X, ' - ')/) 

920  FORMAT( 12X,I5, 5F10. 2) 

INDIVIDUAL  TEST  LOOP 

RSDL=0.0 
WTST=0. 0 
DO  700  L 1 = 1 , NTST 

INITIALIZE  ALL  PARAMETERS 


N0=0 

XX= 1 . OOE+25 
DU1T=0. 0 
XNRME 1=0.0 
XNRMS 1=0.0 
DO  100  11=1,6 
EPSB(II )  =  0. 0 
DEPS(II )=0.  0 
SIGB(II )  =  0. 0 
100  DSIGdI  )  =  0. 0 
U=0.0 
DELU=0.0 
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C  ASSIGN  VALUES  TO  THE  TEST  UNDER  CONSIDERATION 
C 

DO  105  11=1,19 
105  PROP(II )=PRP1 (II ) 

PROP(20)=PRP2(L1, 1 ) 

PR0P(21)=PRP2(L1,2) 

PATH* PROP (8) 

PC0N=PRP2 (L 1 , 3 ) 

KINDsIDAT ( 1 ) 

CNFRsRDAT(l) 

ina=iDAT(2) 

ERMX=RDAT(2) 

LARG=IDAT ( 3 ) 

DO  110  11=1,3 
110  SIGB(II)=PCON 
DO  115  11=1,5 
115  PINC(1,II)=0.0 
PINC(1,  2)=PC0N 
C 

C  ACCOUNT  FOR  CHANGES  IN  THOSE 

C  PARAMETERS  WHOSE  OPTIMAL  VALUES  ARE  BEING  SOUGHT 

C 

IF(JOPT.EQ.O)  GO  TO  125 
DO  120  11= 1 ,NOPT 

120  PROP(IOPT(II ) )=XL(II)+(XU(II )-XL(II) )• 

•  (SIN (ZOPT (II ) /SF (II ) ) )**2 

125  CONTINUE 
C 

C  SCALE  PARAMETERS 

C 

PROP(  3)=PR0P(  3)/3.0/SQRT(3.0) 

PROP(  7)=PR0P(  7)*3-0 
PROP(21)=PROP(21)«3.0 
C 

C  LOADING  HISTORY  SEGMENT  LOOP 

C 

IF( JCOD.EQ. 1 )  WRITE (6, 900)  LI 
DO  600  L2=1,NSEG(L1 ) 

DO  130  11=1,6 
ILOD(II )=LTYP(L1 ,L2, II) 

130  VLOD (II ) =VALU (L 1 , L2, II ) 

NINCsLTYP (L 1 , L2 , 7 ) 

SRATsVALU (L 1 , L2 , 7 ) 

C 

C  DETERMINE  FIRST  INCREMENTS 

C 

DlsFLOAT(NINC) 

D1=1.0/D1 

IP (SRAT . NE . 1 . 0 )  D 1 = ( 1 . 0-SR AT ) / ( 1 . 0-SRAT»«NINC ) 

DO  135  11=1,6 
D2=VL0D(II )-SIGB(II ) 

IF(ILODdl).EQ.I)  D2«VL0D(II)-EPSB(II) 

135  DL0D(II)=D1»D2 
C 
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C  CHANGE  THE  SIGN  OF  THE  STRAIN  ESTIMATE  AT  THE 
C  BEGINNING  OF  THE  NEXT  SEGMENT  IF  THERE  HAS  BEEN 
C  UNSTABLE  BEHAVIOR  AT  THE  END  OF  THE  PREVIOUS  ONE 
C 

DELUs  0. 01*DELU 
DO  140  11=1,6 
DSIG(II )=  0. 01*DSIG(II ) 

140  DEPS (II )=-0. 01*DEPS (II ) 

C 

C  INCREMENT  LOOP 
C 

DO  500  INCR=1,NINC 
C 

C  ITERATION  LOOP  (SUCCESSIVE  APPROXIMATION) 

C 

DO  400  ITNO= 1 , ITMX 
C 

C  ESTABLISH  INCREMENTAL  PROPERTIES 

C 

NJ=N0+1 

K7=ITN0 

CALL  CL AY ( 3 , N J , K7 , PROP , STOR , SIGB , EPSB , 

•  DSIG , DEPS , C , U , DELU , GAMA , KIND , LARG ) 

C 

C  FORM  AND  MODIFY  STIFFNESS  MATRICES 

C 

DO  200  11=1,3 
S(7,II)=GAMA 
S(7,II+3)=0.0 
S(II,7)=1.0 
200  S(II+3,7)s0.0 
S(7,7)=-1.0 
R(7)=0. 0 
DO  210  11*1,6 
DO  205  JJ=1,6 
205  S(II,JJ)=C(II,JJ) 

R(II )=DLOD(II ) 

IF(ILODdl).EQ.O)  GO  TO  210 
S(II,II)*XX 
R(II )=DLOD(II )*XX 
210  CONTINUE 
C 

C  SOLVE  FOR  STRAIN  INCREMENT 
C 

NN=7-KIND 
DO  230  II=1,NN 
D2= 1.0/S (II, II) 

R(II)=R(II)*D2 
DO  215  JJ=II,NN 
215  S(II,JJ)=S(II,JJ)«D2 
IF(II.EQ.NN)  GO  TO  230 
IL*II+1 

DO  225  JJ=IL,NN 
D2=-S(JJ,II) 


* 
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DO  220  KK=II,NN 

220  S ( J J ,  KK )  =S  ( J  J ,  KK )+D2*S ( II , KK ) 

225  R (JJ )=R ( JJ )*D2*R (II ) 

230  CONTINUE 
ICsNN 

DO  235  11=2, NN 

ICsIC-1 

ILsIC+1 

DO  235  JJ=IL,NN 
235  R(IC)=R(IC)-S(IC,JJ)»R(JJ) 

COMPUTE  STRESS  INCREMENT 

DO  245  11=1,6 

D2=0.0 

IF(II.LT.4)  D2=R(7) 

DO  240  JJ=1,6 
240  D2=D2+C(II,JJ)»R(JJ) 

245  RP(II)=D2 

COMPUTE  ERROR  NORMS 

ESIG=0. 0 
EEPS=0.0 
D1=0.0 
D2=0. 0 

DO  250  11=1,6 

ESIG=ESIG+ABS (RP(II )-DSlG(II ) ) 

EEPS=EEPS+ABS (R  (II)-DEPS(II)) 
D1=DWABS(RP(II)) 

250  D2=D2+ABS(R  (II)) 

DU1P=D1 

IF(DUIP.LT.DUIT)  DU1P=DU1T 
ESIG=ESIG/DU 1 P 
EEPS=EEPS/D2 

CHECK  FOR  CONVERGENCE 

IF(ESIG.LT.ERMX.AND.EEPS.LT.ERMX)  GO  TO  405 

APPLY  AITKEN'S  NORM  ACCELERATION 

CNFS =1.0 
CNFE= 1.0 

IF(ITNO.EQ. 1.0R. ((-1 )**ITNO ).GT.O)  GO  TO  255 
CALL  AITKEN (XNRMS2, XNRMS 1 , D 1 , CNFS , CNFR ) 

CALL  AITKEN ( XNRME2 , XNRME 1 , D2 , CNFE , CNFR ) 

255  D1=0. 0 

D2=0.0 

DO  260  11=1,6 

DSIG ( II ) *CNFS »RP ( II )♦ ( 1 . O-CNFS ) 4DSIG ( II ) 
DEPS(II)=CNFE*R  (IlMl.O-CNFE)»DEPS(II) 
D1=D1+ABS(DSIG(II ) ) 

260  D2=D2+ABS (DEPS (II ) ) 
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IF (KIND . EQ . 0 )  DELU=CNFS*R ( 7 )♦ ( 1 . O-CNFS )*DELU 
XNRMS2=XNRMS 1 
XNRME2=XNRME1 
ZNRMS 1 =D 1 
XNRME1=D2 
400  CONTINUE 

IF(JCOD.EQ.I)  WRITE (6, 904)  NO+1 ,ESIG,EEPS 
GO  TO  605 
C 

C  UPDATE  TOTAL  VALUES 
C 

405  NOsNO+1 
DU1T=0.0 
DO  410  11=1,6 
DLOD ( II ) =DLOD (II )  »SRAT 
DSIG(II)=RP(II) 

DEPS(II)=R(II) 

SIGB(II )=SIGB(II )+DSIG(II ) 

EPSB(II )=EPSB(II )+DEPS(II ) 

410  DU1T=DU1T+0. 1*ABS (SIGB(II)) 

IF(KIND.EQ.O)  DELU=R(7) 

U=U+DELU 

C 

C  STORE  INCREMENTAL  VALUES  FOR  FUTURE  PLOTTING 

C 

NJsNO+1 

PINC(NJ , 1 )= (SIGB(3)-SIGB( 1 )) 

PINC (NJ , 2 ) = (SIGB ( 1 )4SIGB ( 2 )*SIGB ( 3 ) )  /3- 0-D 
PINC(NJ,3)=  U 

PINC (N J , 4 ) = ( EPSB ( 1 )+EPSB ( 2 )+EPSB ( 3 ) ) • 1 00. 00 
PINC(NJ, 5)=  EPSB(3)#100.00 
C 

C  PRINT  INCREMENTAL  VALUES 

C 

IF(JCOD.EQ. 1 )  WRITE (6, 908)  NO, (EPSB(II),II=1,6), 

•  (SIGB(II),II=1,6),U,ITN0 

500  CONTINUE 

IF(JCOD.EQ. 1)  WRITE(6,912) 

600  CONTINUE 
605  CONTINUE 

IF(JCOD.EQ.O)  GO  TO  615 
WRITE (6, 9 16)  LI 
DO  610  11=2, NJ 
JJ=II-1 

610  WRITE (6, 920)  JJ,PINC(II,5), (PINC(II,KK),KK=1,4) 

615  CONTINUE 
C 

C  COMPUTE  RESIDUALS  IF  MODEL  CALIBRATION  IS  REQUIRED 
C 

ITST=L1 

IF(JCOD.EQ. 1)  GO  TO  620 
CALL  RSDUAL (NJ , ITST , RTST ) 

RSDL=RSDL+RTST*W1 (LI ) 

WTST=WTST+W 1 (L 1 ) 


620  CONTINUE 


PLOT  MODEL  PREDICTIONS  IP  PLOTTING  IS  REQUIRED 

IF(JCOD.EQ. 1. AND. JPLT.BQ. 1 )  CALL  PLTCHK (NJ , ITST ) 
700  CONTINUE 

IF(UTST.NE.O.O)  RSDL=RSDL/WTST 

RETURN 

END 


SUBROUTINE  AITKBN(X2,X1,XX,XC,XL) 

THIS  SUBROUTINE  COMPUTES 

THE  AITKEN'S  CONVERGENCE  FACTORS. 

XCsl.0 

DX=2.0»I1-X2-XX 
IF(DX.EQ.O.O)  RETURN 
XC=(X  *  C?)/DX 
IF(XC.LT.XL)  XC=XL 
IF(XC.GT. (1.0/XL))  XC»1 . O/XL 
RETURN 
END 


SUBROUTINE  PLTCHK (NP, IT) 

THIS  SUBROUTINE  SERVES  AS  THE  DRIVING  PROGRAM 
FOR  THE  ASSOCIATED  PLOTTING  SUBROUTINES  JSDPLT, 
LSHIFT  AND  BORDER. 

WRITTEN  BY  J.S.  DE  NATALE, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

VERSION  I:  MAY  1982. 

C0»«0N/BLKD/NPLT(6 ) , KPLT (6, 3 ) , NEXP ( 6 , 3 ) , 

*  XV(6,3,20),YV(6,3,20),PINC(61,5) 
DIMENSION  IPLT(IO) 

CHECK  WHETHER  PLOTTING  IS  REQUIRED  (NOTE: 

ADDITIONAL  RELATION  TYPES  MAY  BE  INCLUDED  AS  NEEDED) 

IP=0 

DO  100  11*1,10 
100  IPLT(II):0 

DO  110  II:1,NPLT(IT) 

110  IPLT (KPLT (IT, II ) ): 1 
IF(IPLT(  D.BQ.1) 

•  CALL  JSDPLT (IT, IP, NP,  1,  2,»  Q  V  P  » 
IF (IPLT (  2).EQ. 1 ) 


o  o  o  ooonnoono  ono 


•  CALL  JSDPLT(IT,IP,NP,  1,5,'  Q  ','EPS-Z') 

IF(IPLI(  3).EQ. 1 ) 

•  CALL  JSDPLT (IT, IP, NP,  2,  5,*  P  ','EPS-Z') 

IP(IPLI(  4).EQ. 1) 

•  CALL  JSDPLT (IT, IP, NP,  3,  5, '  0  ','EPS-Z') 

IP(IPLT(  5).EQ. 1 ) 

•  CALL  JSDPLT (IT, IP , HP,  A,  5 , • EPS-V • , • EPS-Z • ) 
RETURN 

END 


SUBROUTINE  JSDPLT (IT , I P , NP , IT , IX , TY , TX ) 

THIS  SUBROUTINE  DIRECTS  THE  PLOTTING  OF  THE  SPECIFIED 
EXPERIMENTAL  OBSERVATIONS  AND  THEORETICAL  MODEL  PREDICTIONS. 

WRITTEN  BY  J.S.  DE  NAT ALE, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

VERSION  I:  MAY  1982. 

COWK)N/BLKD/NPLT  (6 ) ,  KPLT  ( 6 , 3 ) ,  NEXP  (6 , 3 ) , 

'  XV(6,3,20),YV(6,3*20),PINC(61,5) 

DIMENSION  XLAB( 11), YLAB( 1 1 ) , ICPX(61 ) , ICPY (61 ) , ICEX (20) , ICEY (20) 
DIMENSION  FSPA(101),F0RM(4),FSYM(2) 


DOUBLE  PRECISION  TX,TY 


DATA 

FSPA/' 

V  IX. 

','21, 

V  3X, 

'.’  *1, 

.’  5X, ' 

6X,« 

11,' 

.'  8X, 

V  9X, 

• , ’ iox, 

'.MIX, 

'.'12X, 

. ' 13X, ' 

, '  14X, ' 

,’15X,' 

,'16X, 

'.’17X, 

* » ' 18X, 

' » ' 19X, 

' , ' 20X, 

,’21X,' 

,'22  X,’ 

,’23X,’ 

» '24X, 

V25X, 

'.'26X, 

' . '271, 

'.'28X, 

,'29X,’ 

,'30X,' 

,’31X,' 

,'32X, 

’,'331, 

' » '34X, 

V35X, 

'.'361, 

. '37X, ' 

,’38X,' 

,  '39X, ' 

.'40X, 

','411, 

',’*21, 

’,'43X, 

*.'44X, 

,’45x,' 

» '46X, ' 

,’471,' 

,  '48X, 

' , '49X, 

' , '50X, 

'.'Six, 

'.’52X, 

,’53X,' 

.’5AX,’ 

,  '55X, ' 

,’56X, 

' , '57X, 

*,'58x, 

’,'591, 

' . '60X, 

,’61X,’ 

,’62X,’ 

,’631,' 

,’64X, 

' , '65X, 

','661, 

' , '67X, 

'.’68X, 

, '69X, ' 

,’70X,’ 

,’71X,' 

» '72X, 

’,’131, 

’,'7*1, 

' , '75X, 

'.’76X, 

,'Hl,' 

,’78X,» 

,  *791,  ’ 

,'80X, 

'.'BIX, 

' , '82X, 

'.'83X, 

',’84X, 

,’85 X,’ 

,  '86X, ' 

,’87X,' 

,’88X, 

',,89X, 

'.'90X, 

’.'91X, 

’.'92X, 

,’93X,' 

,’94X,’ 

,’95X,’ 

,'96X, 

V97X, 

•,’98x, 

' , '99X, 

’.’99X, 

/ 

DATA 

FORM/' ( 1H+ 

'.'19X, 

«  « 

» 

t  9 
f 

’/ 

DATA 

FSYM/'1H«) 

* , ' 1H#) 

»/ 

IP=IP-*-1 

NE=NEXP(IT, IP) 

ESTABLISH  MINIMUM  AND  MAXIMUM  AXES  VALUES 

XMINspINC ( 1 , IX ) 

XMAX-PINC ( 1 , IX ) 

YMIN=PINC ( 1 , IY ) 

YMAX=PINC( 1 , IT) 

DO  100  IIs2,NP 

IF(XMIN.GT.PINC(II,IX))  XMIN=PINC(II,IX) 
IF(XMAX.LT. PINC(II, IX ) )  XMAX«PINC(II,IX) 
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IF ( THIN . GT . PINC (II , IY ) )  YMIN=PINC ( II , IT ) 

100  IF(YMAX.LT.PINC(IX,IY) )  YMAX*PINC(II,IT) 
IF(NE.BQ.O)  GO  TO  120 
DO  110  11=1, NE 

IF(XMIN. GT. XV (IT, IP, II ) )  XMIN=XV (IT, IP, II ) 
IF(XMAX.LT.XV(IT,IP,II))  XMAX=XV(IT,IP,II) 
IF(TMIN.GT.TV(IT,IP,II))  YHIN=YV(lt,IP,II) 

110  IF(YHAX.LT.YV(IT,1P,II))  TMAX=TV(IT,IP,II) 

120  CONTINUE 

ESTABLISH  ADJUSTED  MAXIMUM  AND  MINIMUM  AXES  VALUES 

CALL  LSHIFT(XMIN.XMAX) 

CALL  LSHIFT (TMIM , TMAX ) 

ESTABLISH  AXES  LABLES 

XDST=XMAX-XMIN 
YDST=YMAX-YMIN 
XINC=XDST/10.0 
YINC=TDST/10. 0 
XLAB(  1  )=XMIN 
XLAB(11 )=XMAX 
YLAB(  1 )=TMIN 
YLAB ( 1 1 ) =YMAX 
DO  130  11=2,10 
AA=FLOAT(II) 

AA=FL0AT(II-1) 

XLAB(II )=XMIN+XINC«AA 
130  TLAB ( II ) =YM1N+YINC*AA 

ESTABLISH  DATA  COORDINATES 

DO  140  11=1, NP 

A1= ( (PINC (II , IX )-  XMIN)/XDST)» 100.0 

IF(AI.LT.O.O)  A1=0.0 

A2=AINT(A1) 

IF(A1-A2.GT.0.5)  A2=A2+1.0 
ICPX(II )=IFIX(A2) 

A 1=( (PINC (II, IT)-  TMIN)/TDST)»  50.0 
IF(AI.LT.O.O)  A1=0. 0 
A2=AINT(A1) 

IF(A1-A2.GT.0.5)  A2*A2+1.0 
140  ICPT(II )=IFIX (A2 ) 

IF(NE.EQ.O)  GO  TO  160 
DO  150  11=1, NE 

A1=((XV(IT,IP,II)-XMIN)/XDST)«100.0 
IF(AI.LT.O.O)  A1=0. 0 
A2=AINT(A1 ) 

IF(A1-A2.GT.0.5)  A2=A2*1.0 
ICEX (II ) =IFIX (A2 ) 

A1=((TV(IT,IP,II)-TMIN)/TDST)»  50.0 
IF(AI.LT.O.O)  A1=0.0 
A2=AINT(A1 ) 
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IFU1-A2.GT.0.5)  A2=A2+1 . 0 
150  ICEY<II)=IFIX(A2) 

PRINT  PLOT 

160  WRITE (6, 170) 

170  FORMAT ( 1H 1 / ) 

11=1 

JJ=0 

DO  210  12=1,51 

PRINT  HEADINGS  AND  AXES  LABLES 
K7=I  2 

CALL  BORDER ( II , J J , K7 , XLAB , YLAB , TX , TY ) 
PRINT  DATA  POINTS 


13=51-12 
DO  180  J2=1,NP 

IF(ICPY(J2).EQ.I3)  FORM(4 )=FSYM( 1 ) 
IF(ICPY(J2).EQ.I3)  P0RM(3)=FSPA(ICPX(J2)+1 ) 
180  IF(ICPY (J2).EQ. 13)  WRITE(6,FORM) 

IF(NE.EQ.O)  GO  TO  200 
DO  190  J2=1,NE 

IF(ICEY(J2).EQ.I3)  F0RM(4 )*FSYM(2) 
IF(ICEY(J2).EQ.I3)  FORM(3)=FSPA(ICEX(J2)+1 ) 
190  IF(ICEY(J2).EQ.I3)  WRITE(6,F0RM) 

200  CONTINUE 
210  CONTINUE 

WRITE(6,220)  (XLAB(MM),MM=1 , 11) 

220  FORMAT(/l6X, 1 1 (F6.4 , 2X) ) 

WRITE (6, 230)  TX 
230  FORMAT ( /20X , 4 8X , A5 ) 

RETURN 

END 


SUBROUTINE  LSHIFT(PMIN, PMAX) 

THIS  SUBROUTINE  ADJUSTS  THE  MINIMUM 
AND  MAXIMUM  AXES  VALUES,  AS  REQUIRED 
BY  SUBROUTINE  JSDPLT. 

DX=(PMAX-PMIN)/10.0 

XL=PMIN-DX 

XU=PMAX+DX 

DO=ABS(PMIN) 

D1=ABS(PMAX) 

IF(DI.LT.DO)  D1=D0 
NN=21 

DO  100  JJ=1 ,40 
NN=NN-1 


o  o  o  o  ooo 


C1*10.0**NN 
IF(DI.GE.CI)  GO  TO  110 
100  CONTINUE 
110  DUPHAX/C1 
II=IFII(D1 ) 
DU=FL0AT(II)«C1 
IF(DU.LT.PMAX)  11*11+1 
DU=FL0AT(II )*C1 
IF(DU.LE.XO)  GO  TO  120 
C1=C1/10.0 
GO  TO  110 
120  D1sPMIN/C1 
II=IFIX(D1) 
DL=FLOAT(II)«C1 
IF(DL.GT.PMIN)  11=11-1 
DL=FLOAT(II)»C1 
IF(DL.GE.XL)  GO  TO  130 
C1=C1/10.0 
GO  TO  110 
130  PMINrDL 
PHAX=DU 
RETURN 
END 


SUBROUTINE  BORDER (II , J J , 12, XLAB, TLAB, XT , YT ) 

THIS  SUBROUTINE  PRODUCES  THE  PLOT  HEADINGS 

AND  AXES  LABLBS ,  AS  REQUIRED  BT  SUBROUTINE  JSDPLT. 

DIMENSION  XLAB ( 1 1 ) t  TLAB (11) 

DOUBLE  PRECISION  XTfYT 
K=0 

IF(I2.EQ. 1.0R.I2-II.EQ.5)  K=1 
IF(KK.EQ. 1 )  11*12 

IFUOC.BQ.  1)  JJsJJ+1 
LL*12-JJ 

IF(I2.EQ.  1)  WRITE(6, 110)  TLAB(LL) 

IF(I2.EQ.  1)  GO  TO  100 
IF(I2.EQ.26)  WRITE (6, 120)  YT, YLAB(LL) 

IF(I2.EQ. 26)  GO  TO  100 
IF(I2.EQ.51)  WRITE(6, 110)  YLAB(LL) 

IF(I2.EQ.51 )  GO  TO  100 
IF(RX.EQ.  1)  WRITE(6,130)  YLAB(LL) 

IFOCX.EQ.  0)  WRITE (6, 140) 

100  CONTINUE 

110  F0RMAT(1X,5X,2X,F10.5,2X,  'I* ,  10(* - 1')) 

120  FORMAT( 1X,A5, 2X,F10.5,2X, ,+-1 ,971, ,-+’ ) 

130  FORMAT(1X,5X,2X,F10.5,2X, »«~»,97X, »-+' ) 

140  FORMAT ( IX, 5X, 2X,  10X,2X,»I  »,97X,'  I') 

RETURN 

END 
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SUBROUTINE  RSDU AL (NP , IT , RTST ) 


THIS  SUBROUTINE  COMPUTES  THE  WEIGHTED 
ABSOLUTE  OR  SQUARED  RESIDUAL  BETWEEN  THE 
EXPERIMENTAL  OBSERVATIONS  AND  THE  CONSTITUTIVE 
MODEL  PREDICTIONS. 

WRITTEN  BY  J.S.  DE  NAT ALE, 

DEPARTMENT  OP  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

VERSION  I:  MAY  1982 

COMMON/BLICB/W 1  (6  ) ,  W2  ( 6 , 3  ) ,  W3  ( 6 , 3 , 20) 
COMMON/BLKD/NPLT ( 6 ) , KPLT ( 6 , 3 ) , NEXP ( 6 , 3 ) , 

•  XV(6,3,20),YV(6,3,20),PINC(61,5) 
COMMON /BLKE /KRSL , POWR , XMXR ( 5 ) , YMXR ( 5  ) , XMXP (6,3) 

•  YMXP(6, 3),PLIM,TLIM 

INDIVIDUAL  PLOT  LOOP 

RTSTsO.O 

WPLTsO.O 

DO  400  IP=1 ,NPLT(IT) 

IDENTIFY  RELEVANT  MAXIMA 

GOTO  (100,102,104,106,108),  KPLT(IT,IP) 

100  IX=2 
IY=1 

XALLsXMXR(l) 

YALLsYMXR ( 1 ) 

GO  TO  120 
102  IX=5 
IY=  1 

XALLsXMXR ( 2 ) 

YALLsYMXR (2) 

GO  TO  120 
104  IX=5 
IY*2 

XALL=XMXR(3) 

YALLsYMXR (3) 

GO  TO  120 
106  IX=5 
IT=3 

XALLsXMXR (4) 

YALLsYMXR (4) 

GO  TO  120 
108  IXs5 
IYs4 

XALLsXMXR (5) 

YALLsYMXR (5) 

120  CONTINUE 


ooo  n  n  n  non  non  non  non 


COMPUTE  AND  SUM  WEIGHTED  RESIDUALS 

IF(KRSL.NE.O)  LL=1 
IF(KRSL.EQ.O)  LL=2 
RPLTsO.O 
WPNT=0. 0 
XMAX=XMXP(IT,IP) 

YMAX=YMXP(IT,IP) 

CHECK  FOR  ABSOLUTE  OR  RELATIVE  TEST  SCALING 


XMXTsXMAX 

YMXT=IMAX 

IF( (TMAX/YALL).LT.TLIM)  TMZT>TALL*TLIM 

INDIVIDUAL  EXPERIMENTAL  OBSERVATION  LOOP 

DO  300  IE=1,NEXP(IT,IP) 

RPNT=0. 0 

MODEL  PREDICTIONS  LOOP  (IN  CALLED  SUBROUTINE) 


IMIN=0 

IPLT=IP 

IEXP=IE 

IF(KRSL.NE.O)  CALL  VRTICL(n,IPLT,IEXP,NP,lX,IY, 

•  LL,YE,DR,IMIN) 
IF(KRSL.EQ.O)  CALL  EUCLID(IT,IPLT,IEXP,NP,IX,IY, 

•  LL,YE,DR,IMIN,XMXT,YMXT) 

CHECK  FOR  ABSOLUTE  OR  RELATIVE  POINT  SCALING 

IF(IMIN.EQ.O)  GO  TO  200 
YEsABS(YE) 

YMXV=TE 

IF( (YE/YMAX ) .LT. PLIM)  YMXV=YMAX»PLIM 
RPNT* ( (ABS (DR ) /YMXV )• (YMAX/YMXT ) )»»P0WR 
RPLTsRPLT+RPNT»N3 (IT, IP, IE ) 

WPNT=WPNT+N3 (IT , IP, IE ) 

200  CONTINUE 
300  CONTINUE 

RPLT=RPLTAfPNT 
RTST*RTST+RPLT»W2 (IT, IP ) 

WPLT=NPLT+W2(IT,IP) 

400  CONTINUE 

RTST«RTST/WPLT 

RETURN 

END 


SUBROUTINE  EUCLID (IT , IP , IE , NP , IX , IY , LL , YE, DR , 
«  IMIN , XMXT , YMXT ) 
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THIS  SUBROUTINE  COMPUTES  THE  MINIMUM  EUCLIDEAN 
DISTANCE  BETWEEN  THE  EXPERIMENTAL  OBSERVATION  AND 
THE  MODEL  PREDICTION  ■CURVE". 

WRITTEN  BY  J.S.  DE  NAT ALE, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

VERSION  I:  MAY  1982. 

COMMON/BLKD/NPLT ( 6 ) , KPLT ( 6 , 3 ) , NEXP ( 6 , 3 ) , 

•  XV(6,3,20),YV(6,3,20), PINC (61,5) 

DETERMINE  THE  PREDICTION  POINT  WHICH  IS  NEAREST 
TO  THE  EXPERIMENTAL  OBSERVATION  IN  A  SUITABLY  SCALED 
X-Y  SPACE 

ML=0 

DR= 1 . 00E+02 

XE=XV (IT, IP,IE)/XMXT 

YE=YV(IT,IP,IE)/YMXT 

DO  100  KL=LL,NP-1 

XP=PINC(KL,IX)/XMXT 

YPsPINC (KL , IY ) /YMXT 

DV=SQRT ( ( XE-XP ) ••2+ ( YE-YP )««2) 

IF(DR.GT.DV)  MLsKL 
100  IF(DR.GT.DV)  DR=DV 

IF(ML.EQ.  0)  GO  TO  120 

SEARCH  ADJACENT  SEGMENTS  FOR  AN  EVEN  SMALLER  DISTANCE 
NL=ML 

DO  110  MM= 1 , 2 

K1=NL-1 

K2=NL 

X 1 =PINC (X 1 , IX ) /XMXT 

Y1=PINC(K1,IY)/YMXT 

X2=PINC (K2, IX ) /XMXT 

Y2=PINC(K2,IY)/YMXT 

DX=X2-X1 

DY=Y2-Y1 

BETAs  ( (XE-X 1 )  •DX+CYE-T 1  )*DY )  /  (DX”2+DY”2 ) 

IF(BETA.LT.O.O)  BETA=0.0 

IF(BETA.GT.I.O)  BETA=1.0 

XM=XUBETA»DX 

YM=YUBETA»DY 

DVsSQRT( (XE-XM)#t2+(YE-IM)**2) 

IF(DR.GT.DV)  DR=DV 
110  NL=NL+1 
LL=ML 
IMINr 1 

120  DR=DR*YMXT 
YE = YE* YMXT 
RETURN 
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END 


SUBROUTINE  VRTICL ( IT , IP , IE , NP , IX , I Y , LL , YE , DR , IMIN ) 

THIS  SUBROUTINE  COMPUTES  THE  MINIMUM  VERTICAL 
DISTANCE  BETWEEN  THE  EXPERIMENTAL  OBSERVATION  AND  THE 
MODEL  PREDICTION  "CURVE*. 

WRITTEN  BY  J.S.  DE  NAT ALE, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

VERSION  I:  MAY  1982. 

COMSON/BLKD/NPLT ( 6 ) , KPLT (6,3),NEXP(6,3), 

•  XV(6,3,20),YV(6,3,20),PINC(61,5) 

DETERMINE  THE  PREDICTION  POINTS  WHICH  MOST  CLOSELY 
BOUND  (WRT  THE  X-VALUE)  THE  EXPERIMENTAL  OBSERVATION 

XE=IV(IT,IP,IB) 

YE=YV (IT, IP, IE ) 

DO  100  EL*LL,NP-1 
KU=KL+1 

XLsPINC (XL , IX ) 

YL=PINC(KL,IY) 

XU=PINC(KU,IX) 

YU=PINC(KU,IY) 

IF(XE.EQ.XL)  GO  TO  110 
IF(XE.EQ.XU)  GO  TO  120 
IF(XE.LT.XL. AND.XE.GT.XU)  GO  TO  130 
100  IF(XE. GT.XL.AND.XE.LT. XU)  GO  TO  130 
GO  TO  150 
110  YPsYL 

GO  TO  140 
120  YPrYU 

GO  TO  140 

130  BETA= ( XE-XL ) / ( XU -XL ) 

YP=TL-»-  BETA"  (YU -TL ) 

140  DR=ABS(YE-YP) 

LL=XL 
IMINs 1 

150  CONTINUE 
RETURN 
END 


SUBROUTINE  BSMOPT(XX) 

THIS  SUBROUTINE  DIRECTS  THE  CALIBRATION  OF  THE 
BOUNDING  SURFACE  PLASTICITY  MODEL.  THE  OPTIMAL 
VALUES  OF  THE  MATERIAL  AND  MODEL  PARAMETERS  ARE 
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SELECTED  BY  MINIMIZING  THE  WEIGHTED  RESIDUAL  BETWEEN 
THE  EXPERIMENTAL  OBSERVATIONS  AND  THE  MODEL  PREDICTIONS 

WRITTEN  BY  J.S.  DE  NAT ALE, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

VERSION  I:  MAY  1982. 

COMMON /BLK1/ETAF, ETAX , DMIN , FDER 
COMMON /BLK2/NDIM, NFMX , C 1 , C2, C3 
COMMON/BLK3/NI,NF,NG,Fl,F2,S1,S2,X1(17),X2(17) 

COMMON /BLK4 /EPFA , EPS A , ICON , JCON , I PRN , ICUB , XE ( 1 7 ) 
COMMON/BLK5/D1(17),D2(17),G1(17),G2(17),DG(17),DX(17), 

»  HI (17, 17),H2<17, 17) 

COMMON/BLK6/EPFL , EPXL , NFLX 
COMMON/BLK7/SFUN , SF ( 1 7 ) , XL ( 1 7 ) , XU ( 1 7 ) 

DIMENSION  XX (17) 

C  READ  CALIBRATION  SPECIFICATIONS 

CALL  CALDAT 

INITIALIZE  GLOBAL  SEARCH 

DO  100  II=1,NDIM 
TX=X2(II) 

TD=X2 ( II ) +DX ( II ) 

X2(II )=ASIN(SQRT( (TX-XL (II ) )/(XU(II )-XL(II ) ) ) ) 

DX(II )=ASIN (SQRT( (TD-XL(II ) )/(XD (II)-XL(II) ) ) )-X2(II) 
X2(II )=X2(II )*SF(II ) 

100  DX(II)=DX(II)«SF(II) 

NF=0 
NGsO 
NI=0 

F2=FFUN(X2) 

NF=NF+1 

DO  110  II=1,NDIM 
IG=II 

110  G2(IG)=GFUN(IG,0,NDIM) 

NG=NG+1 

CALL  MVPR0D(H2,G2,D2,NDIM) 

S2=DTPR0D(G2,D2,NDIM) 

IF(IPRN.GE.I)  CALL  PRNOUT(I.NDIM) 

PERFORM  LINE  SEARCH 

200  CALL  SEARCH 

IF  NO  MINIMUM  HAS  BEEN  FOUND,  UPDATE  THE 
INVERSE  HESSIAN,  COMPUTE  A  NEW  SEARCH  DIRECTION. 

AND  PROCEED  WITH  THE  NEXT  ITERATION 

IF(JCON.EQ.I)  GO  TO  300 
CALL  BFGSUP (S2 , IHES , NI , NDIM ) 
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JCONrO 

IF(ICON.EQ. 3)  CALL  CONCUR ( DF , S2 , DX , NDIM ) 
IF(JCON.EQ.I)  CALL  PRNOUT( 1.NDIM) 
IF(JCON.EQ.I)  CALL  PRN0UT(5,NDIM) 
IF(JCON.EQ.I)  GO  TO  300 
GO  TO  200 

TERMINATE  GLOBAL  SEARCH 

300  CONTINUE 

DO  310  11=1, NDIM 
310  XX(II)=X2(II) 

RETURN 

END 


SUBROUTINE  CALDAT 

THIS  SUBROUTINE  READS  IN  ALL  ADDITIONAL 
INFORMATION  REQUIRED  TO  DIRECT  THE  MODEL 
CALIBRATION  ALGORITHM. 

WRITTEN  BY  J.S.  DE  NAT ALE, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

VERSION  I:  MAY  1982. 

COMMON /BLKA/NTST , JOPT , JPLT, NOPT , KOPT (IT) 

COMMON/BLKC/PRP 1(19), PRP2 (6,4), IDAT ( 3 ) , RDAT ( 2 ) , 

•  NSEG(6),LTYP(6,4,7),VALU(6,4,7) 

COMMON /BLR  1 /ETAF, ETAX , DMIN , FDER 
C0MM0N/BLK2/NDIH , NFMX , C 1 , C2, C  3 
COMMON/BLK3/NI,NF,NG,F1,F2,S1,S2,X1(17),X2(17) 

COMMON /BLK4 /EPFA , EPS A , ICON , JCON , I PRN , ICUB , XE ( 1 7 ) 
COMMON/BLK5/D1(17),D2(17),G1(17),G2(17),DG(17),DX(17), 

•  H1(17,17),H2(17,17) 

COMMON /BLK6/EPFL , EPXL , NFL I 

COMMON /BLK7/SFUN , SF ( 1 7 ) , XL ( 1 7 ) , XU ( 1 7 ) 

FORMAT  STATEMENTS 

800  FORMAT (1615) 

804  F0RMAT(7E10. 3,215) 

808  FORMAT(8E50. 3) 

812  F0RMAT(4 (2I5,E10. 3) ) 

816  FORMAT (5E 10.3, 615 ) 

900  F0RMATOH1) 

904  F0RMAT(5X, 'INPUT  DATA: '/5X, 'CALIBRATION  CONTROL  CODES:'/ 

•  5X, ' - ’// 

•  15X, 'NUMBER  OF  DIMENSIONS  >',14/ 

•  15X, 'MAXIMUM  NUMBER  OF  FUNCTION  EVALUATIONS  »',I4/ 

•  15X, 'OUTPUT  DATA  PRINT  CODE  =',14/ 

•  15X, 'CONVERGENCE  CRITERION  CODE  >',14/ 


118 
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908 

912 

916 

920 


924 


928 


932 

936 

940 


944 


15X, 'CUBIC  INTERPOLATION  OPTION  CODE 

15X, 'HESSIAN  CONDITIONING  OPTION  CODE 

15X, '#  OF  NON-ZERO  HESSIAN  VALUES  TO  BE  READ 

15X, 'PARAMETER  SCALING  OPTION  CODE 

15X,'#  OF  IMPROVED  STARTING  ESTIMATES  REQUIRED 

15X, 'ABSOLUTE  CONVERGENCE  CRITERION  FOR  F(X) 

15X, 'ABSOLUTE  CONVERGENCE  CRITERION  FOR  S(X) 

1 51 VALUE  OF  THE  LINE-SEARCH  PARAMETER,  SlOtt 
15X, 'VALUE  OF  THE  BRACKET-CHECK  PARAMETER,  ROA 
151, 'VALUE  OF  THE  BRACKET -CHECK  PARAMETER,  TAU 
15X, 'VALUE  OF  THE  LINE  SEARCH  EXIT  CRITERION 
15X,' VALUE  OF  THE  FAIL-SAFE  GLOBAL  CRITERION 
15X, '#  OF  FN-EVALUATIONS  BEFORE  FAIL-SAFE  EXIT 
15X, 'ERROR  BOUND:  MACHINE  REPRESENTATION  OF  XX 
15X, 'ERROR  BOUND:  MACHINE  REPRESENTATION  OF  FX 
15X, 'MINIMUM  VALUE  OF  THE  DIFFERENCE  INTERVALS 
15X, 'MAXIMUM  ERROR  IN  FWD-DIFFERENCE  CALCULATN 
15X, 'SCALING  FACTOR  FOR  THE  OBJECTIVE  FUNCTION 


.1*/ 

,1*/ 

,IV 

,1*/ 

,14/ 

, 1PE10. 3/ 

,  E10.3/ 

,  E10.3/ 

,  B10.3/ 

,  E10.3/ 

,  E10.3/ 

,  E10.3/ 
,14/ 

,  E10.3/ 

,  E10.3/ 

,  E10.3/ 

,  BIO. 3/ 

,  E10.3/) 


FORMAT (/5X, 'INPUT  DATA: '/5X, 'INITIAL  VALUES  OF  THE  FUNCTION  ', 
•VARIABLES  XI, X2, . . . ,XN: '/5X,54('-» )/10X, 1P12E10.2) 

FORMAT (/5X, 'INPUT  DATA: */5X, 'LOWER  AND  UPPER  BOUNDS  ON  THE  ', 
•FUNCTION  VARIABLES  XI ,X2, . . . ,XN: • /5X,62( •-' ) ) 

FORMAT ( 1 OX , 1P12E10.2) 

FORMAT (/5X, 'INPUT  DATA: ’/5X, 'INITIAL  VALUES  OF  THE  ', 

•DIFFERENCING  INTERVALS  D1 ,D2, . . . ,DN: '/5X,58('-' ) 

/10X, 1P12E10.2) 

FORMAT (/5X, 'INPUT  DATA: */5X, 'ABSOLUTE  CONVERGENCE  CRITERIA  FOR  ', 

•THE  FUNCTION  VARIABLES  XI, 12 . XN: '/5X,70('-' ) 

/10X, 1P12E10.2) 

FORMAT (/5X, 'INPUT  DATA: '/5X, 'SCALING  FACTORS  FOR  THE  ', 

•FUNCTION  VARIABLES  XI ,X2, . . . ,XN: '/5X,56( •-• ) 

/10X, 1P12E10.2) 

FORMAT (/5X, 'INPUT  DATA: '/5X, 'INITIAL  VALUE  OF  THE  INVERSE  ', 
•HESSIAN  MATRIX :'/5X, 44 (•-•)) 

FORMAT ( 1 OX , 1P12E10.2) 

F0RMAT(1H1,4X, 'OUTPUT  DATA :' /5X, 'THE  RESULTS  OF  THE  OPTIMIZATION  ', 
•ROUTINE  FOLLOW: »/5X, 47 ('-' )//5X, '#-1' , IX, '#-F' , IX, »#-G* , 

9X, 'FI • ,9X, 'F2' ,9X, »S1 • ,9X, 'S2' ,9X, 'XI ' ,9X, »X2- ,9X, 'X3' , 
9Xt'X4\9X,'X5',9X,'X6'/ 

4X,  3(  IX, ' - ' ),  10(5X, ' - ')/) 

F0RMAT(1H 1,41, 'OUTPUT  DATA: '/5X, 'THE  RESULTS  OF  THE  OPTIMIZATION  ', 
' ROUTINE  FOLLOW : ' /5X , 4  7 ( ' - ' ) //5X , » #-I ' , IX , • #-F ' , IX , ' #-G • , 

9X, 'F1',9X,'F2',9X,'S1',9X,'S2',6X,'X1/X7',6X,'X2/X8', 

6X,  'X3/X9',5X,  'X4/X10'  ,5X,  'X5/X11 ' ,5X,  'X6/X12'/ 

4X,  3(  IX, ' - ' ),  10(5X, ' - ')/) 


READ  CALIBRATION  SPECIFICATIONS 

READ (5, 800)  NDIM,NFMX,IPRN,ICON, 

»  I CUB , IHES , NHES , ISCA , NSLN 

READ (5, 804 )  EPFA.EPSA,  C1,C2,C3,EPXL,BPFL,NFLX 
READ (5, 808)  ETAF,ETAX,DMIN,FDER,SFUN 
READ(5,808)  (X2(II),IIs1,NDIM) 

READ(5,808)  (XL(II),II«1,NDIM) 
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BEAD (5 i 808)  (XU (II), 11=1, NDIM) 

BEAD (5, 808)  (DZ (II), 11= 1 ,HDIM) 

IF(ISCA.EQ. 1 )  RE AD (5, 808)  (SF(II ),II=1 ,NDIM) 
IF(ICON.EQ.I)  READ (5, 808)  (XE(II),II«1,NDIM) 

INITIALIZE  INVERSE  HESSIAN 

DO  100  11=1 , NDIM 
DO  100  JJx1,NDIM 
H2(II,JJ)=0.0 

100  IF(II.EQ.JJ)  H2(II, JJ )=-1 . 0 
IF(NHES.EQ.O)  GO  TO  130 
A1=FL0AT(NHES)/4.0 
A2=AINT(A1 ) 

IFU1-A2.GT.0.0)  A2=A2+1. 0 
NH=IFIX(A2) 

NX=NHES-4*(NH-1 ) 

DO  110  11=1, NH 

IF(II.LT.NH)  READ(5, 812)  ((IR,IC,H2(IR,IC)),JJ=1,  4) 
110  IF(II.EQ.NH)  READ(5|612)  ((IR,IC,H2(IR,IC)),JJ=1,NX) 
DO  120  11=1, NDIM 
DO  120  JJ=1,NDIM 

IF(H2(II, JJ ).NE. 0. 0)  H2(II,JJ)=-H2(II,JJ) 

120  IF (H2 (II ,  J  J  ) .  NE .  0. 0  )  H2(JJ,II)=  H2(II,JJ) 

130  CONTINUE 

COMPOTE  IMPROVED  INITIAL  ESTIMATES 

IF(NSLN.EQ.O)  GO  TO  200 
DO  190  II=1,NSLN 

READ(5,816)  PLAM , PKAP , PCON , PFAL , RI , KR 

IF(RI.EQ.O.O)  RI=2.50 

CALL  NEVTON ( PLAM  , PK AP , PCON , PFAL ,  RI ) 

GO  TO  (140,160),  KR 
140  DO  150  J  J= 1 , NDIM 
150  IF(KOPT(JJ).EQ.  9)  X2(JJ)=RI 
GO  TO  180 
160  RC=0.0 

DO  170  JJ=1,NDIM 
IF(KOPT(JJ).EQ. 12)  KK=JJ 
170  IF(KOPT(JJ).EQ.  9)  RC=X2(JJ) 

IF(RC.EQ.O.O)  RCsPRPI (9) 

X2(KK)sRI/RC 
180  CONTINUE 
190  CONTINUE 
200  CONTINUE 

ASSIGN  DEFAULT  VALUES 

IF(CI.EQ.O.O)  CU9.00E-01 
IF(C2.EQ. 0.0)  C2= 1 . 00E-04 
IF(C3.EQ.0.0)  C3«1.00E-01 
IF(EPXL.EQ.O.O)  EPXL* 1 . O0E-04 
IF(ETAF.EQ.O.O)  ETAF=1 .  00E-04 
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IF(ETAZ.EQ.O.O)  ETAI= 1 . 00E-08 
IF(DMIN.EQ.O.O)  DMIN=5. 00E-03 
IF(FDER.EQ.O.O)  FDER= 1 . 00E-04 
IF(SFUN.EQ.O.O)  SFUN= 1 . 00E+00 
SMAL= 1 . 00E-06 
DO  210  11=1, NDIM 
XT=X2(II)+DX(II) 

IF(XT.LE.XL(II) )  DX(II )=XL(II )-X2(II )+SMAL 
210  IF(XT.GE.XU(II))  DXdI)=XU(II)-X2(II)-SMAL 

PRINT  INPUT  DATA 

VRITE(6f 900) 

WRITE (6 , 904 )  NDIM , NFMX , IPRN , ICON , ICUB , IRES , NHES , 
•  ISCA , NSLN , EPFA , EPSA ,  C1,C2,C3,EPXL, 

■  EPFL , NFLX i ETAX , ETAF , DHIN , FDER , SPUN 

WRITE (6 ,908)  (X2 (II), 11=1, NDIM) 

WRITE(6,912) 

WRITE(6,916)  (XL(II),II=1 ,NDIM) 

WRITE(6,916)  (XU(II),II=1,NDIH) 

WRITE (6, 920)  (DX(II), 11*1, NDIM) 

IPdCON. EQ.  1 )  WRITE (6, 924)  (XE(II),II»1,NDIH) 

IF (ISCA. NE. 0)  WRITE (6, 928)  (SF(II), 11=1, NDIM) 
IF(NHES.EQ.O)  GO  TO  230 
WRITE (6, 932) 

DO  220  11=1, NDIM 

220  WRITE<6,936)  (-H2(II, JJ),JJ*1,NDIM) 

230  IFCNDIM.LE.6)  WRITE(6,940) 

IF(NDIM.GT. 6)  WRITE(6,944) 

RETURN 

END 


SUBROUTINE  NEWTON (PRP 1 , PR P2 , PCON , PFAL , RI ) 

THIS  SUBROUTINE  COMPUTES  AN  INITIAL  ESTIMATE 
OF  THE  SURFACE  SHAPE  PARAMETER  R  BY  MEANS  OF  THE 
MODIFIED-NEWTON  ROOT  FINDING  SCHEME. 

WRITTEN  BY  J.S.  DE  NAT ALE, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

VERSION  I:  MAT  1982. 

NI=0 

RO=RI 

DR=0. 0 

ER= 1 . 0E-06 

ALPH = PFAL /PCON 

BETA=PRP1/(PRP2-PRP1 ) 

A1=ALPH»»BETA 

A2=ALPH««(2.0»BETA) 

100  R1=RI 
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HIsHI+1 

IF(NI.GT.IOO)  RI=R0 
IF(HI.GT.IOO)  RETURN 

1 1 0  FF*2. 0» (A 1-A2J+R 1  •  ( A2-2 . 0 )+2 . 0»R 1  M2-R  1  ««3 
GG=A2-2. 0+4 . 0*R 1 -3. 0»R 1 »»2 
HH=4.0-6.0»R1 
UFsFF/GG 

UG=1.0-FF*HH/GG«2 

DRsUF/UG 

RUR1-DR 

IF(ABS(DR).LE.ER)  GO  TO  120 
GO  TO  110 

IF  THE  COMPUTED  VALUE  OF  R  IS 
BETWEEN  1.50  AND  3.50,  THEN  STOP 

120  IF(R1.GE. 1.5)  GO  TO  130 
RIsRI+1.0 
GO  TO  100 

130  IF(R1.LE.3.5)  GO  TO  140 
RI=RI-0. 5 
GO  TO  100 
140  RI=R1 
RETURN 
END 


SUBROUTINE  SEARCH 

THIS  SUBROUTINE  DIRECTS  THE  LINE-SEARCH 
COMPONENT  OF  THE  GLOBAL  MINIMIZATION  ALGORITHM. 

THE  LOGIC  IS  PATTERNED  AFTER  THAT  OF  FLETCHER  (I960). 

WRITTEN  BY  J.S.  DE  NATALE, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

VERSION  I:  MAY  1982. 

C0MM0N/BLK2/NDIM,NFMX,C1,C2,C3 

COMMON/BLR3/NI,NF,NG,F1,F2,S1,S2,X1(17),X2(17) 

COMMON /BLK4 /EPF A , EPSA , ICON , JCON , IPRN , ICUB , XE ( 1 7 ) 
COMMON/BLK5/D1(17),D2(17),G1(17),G2(17),DG(17),DX(17), 
•  H1(17,17),H2(17,17) 

COMMON /BLK6/EPFL , EPXL , NFLX 
COMMON/BLK7/SFUN,SF(17),XL(17),XU(17) 

DIMENSION  DM(17),XM(17) 

INITIALIZE  AND  UPDATE  PARAMETERS 

FEND&F2 
NIsNI+1 
AlxO.O 
A2« 1 . 00E+06 
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F1=F2 

S1=S2 

DO  100  11= 1 ,NDIM 
D1  (II  )=D2(II) 

G1  (II  )sG2(II ) 

X1(I1)=X2(II) 

DO  100  JJ=1,NDIM 
100  H1(II,JJ)=H2(II,JJ) 

ESTABLISH  AN  INITIAL  VALUE  OF  ALPHA 

IF(NI.EQ.  1)  DF=F1 

Z1=10.0«EPFA 

IF(DF.LT.ZI)  DFsZI 

AA=  1.0 

Z1s-2.0«DF/S1 

IF(AA.GT.ZI)  AA=Z1 

IF (LEND . EQ . 1 . AND . LP AS . EQ .  0 )  AA= 1 . 0 

EVALUATE  THE  FUNCTION  F(X1,X2, . . . ,XN), 

BUT  TERMINATE  THE  LINE  SEARCH  IF  THE  CHANGE  IN 
X1,X2, . . . ,XN  IS  SMALL 


LPAS=0 

200  IF(NF.GE.NFMX)  GO  TO  500 
LEND= 1 

DO  210  IIsl.NDIM 

X3*XL  (II  )«■  (XU  (II  )-XL  (II )  )»(SIN  (X2(II )  )  )”2 
DX(II)=AA*D1 (II) 

X2(II )=X1 (II )+DX(II ) 

I4=XL(IIMXU(II)-XL(II))»(SIN(X2(II)))«*»2 

XC=ABS((X3-X4)/X3) 

210  IF(XC.GT.EPXL)  LEND=0 

IF(LEND.EQ. 1.AND.LPAS.EQ. 1)  GO  TO  420 
F2=FFUN(X2) 

NF=NF+1 


CHECK  THE  FAIL-SAFE  EXIT  CRITERION 
AND  TERMINATE  THE  GLOBAL  SEARCH  IF  THE 
FUNCTION  IS  NOT  BEING  SIGNIFICANTLY  REDUCED 

IF( (F1-F2).GT.EPFL)  IX=NF 
IF( (NF-IX).GT.NFLX)  GO  TO  500 
IF(LEND.EQ.I)  GO  TO  300 

CHECK  TO  SEE  THAT  ALPHA  IS  WITHIN  THE 
ACCEPTABLE  INTERVAL.  CHECK  UPPER  BOUND 

DF*F1-F2 

Z1*  -C2*AA*S 1 

IF(DF.GE.ZI)  GO  TO  300 

CHECK  FAILS.  COMPUTE  AN  ACCEPTABLE  VALUE  BY 
QUADRATIC  INTERPOLATION 
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IFCIPRN.EQ.3)  CALL  PRNOUTC 1 ,NDIM) 

ZUAA-A1 

Z2«2 . 0« ( 1 . O+DP/ (Z 1 «S 1 ) ) 

AH*AWZ1/Z2 

C 

C  UPDATE  BRACKET 

C  ENSURE  THAT  ALPHA  IS  NOT  NEAR  THE  ElTRBffiS 
C 

A2rAA 

Z2*C3»Z1 

ir(AH.LT.  (AUZ2))  AH«AUZ2 
IFCAH.GT. (A2-Z2) )  AH*A2-Z2 
AA=AH 
GO  TO  200 
C 

C  CHECK  SUCCEEDS.  EVALUATE  THE  GRADIENT  AND  SLOPE 

C 

300  DO  310  Ild.NDIM 
IG=II 

DMCII )>DZ(II ) 

XM(II)*I2(II) 

310  G2(IG)sGFUN(IG,0,NDIM) 

NGrNG+1 

S2«DTPR0D (G2 , D 1 , NDIH ) 

FH=F2 

SM«S2 

IF(LEND.EQ.I)  GO  TO  400 
C 

C  CHECK  TO  SEE  THAT  ALPHA  IS  VITHIN  THE 
C  ACCEPTABLE  INTERVAL.  CHECK  LOWER  BOUND 
C 

L PAS si 

IFCS2.GE. (C 1 *S 1 ) )  GO  TO  350 
C 

C  CHECK  FAILS.  COMPUTE  AN  ACCEPTABLE  VALUE 

C  BT  QUADRATIC  EXTRAPOLATION 

C 

IFCIPRN.EQ.3)  CALL  PRNOUT ( 1 , NDIM ) 

Z 1 sAA-A 1 
Z2=S2/(S1-S2) 

AHsAAfZ1«Z2 

C 

C  UPDATE  BRACKET 

C  ENSURE  THAT  ALPHA  IS  NOT  NEAR  THE  EXTREMES 

C 

Z2*C3*Z1 

IFCAH.LT. CAA+Z2))  AH«AA+Z2 

z2>9.o«zi 

IFCAH.GT. CAA+Z2))  AH*AA+Z2 

Z2«(A2-AA)/2.0 

IFCAH.GT. (AA+Z2))  AH«AA+Z2 

AUAA 

AAsAH 


o  o  o  noon  o  o  o  o  o  o  o  o 


F1=F2 
S1*S2 
GO  TO  200 

RE-CHECK  LOWER  BOUND  IF  THE 

STRICT  TERMINATION  CRITERION  IS  EMPLOYED 

350  IF(ICDB.EQ.I)  GO  TO  400 

IF(ABS(S2).LE. (-C1*S1)}  GO  TO  400 

CHECK  FAILS.  COMPUTE  AN  ACCEPTABLE  VALUE 
BY  CUBIC  INTERPOLATION 

IF(IPRN.EQ.3)  CALL  PRNOUT(I.NDIM) 

ZUAA-A1 

DD= (2. 0*(F1-F2)+Z1*(S1+S2) )/(Z1**3) 
CCs(S2-SU3.0»DD»(A1«»2-AA»»2))/(2.0»Z1) 

BBs  (  AA«S  1  -A  1  «S2i-3 . 0«DD»  A 1  *AA«Z  1 )  /Z 1 
AH=BB/ ( -CC-SQRT (CC*CC-3* 0*BB*DD) ) 

UPDATE  BRACKET 

ENSURE  THAT  ALPHA  IS  NOT  NEAR  THE  EXTREMES 

A2=AA 

Z2=C3»Z1 

IF(AH.LT.  (AWZ2))  AH*A1+Z2 
IF(AH.GT. (A2-Z2))  AH=A2-Z2 
AA=AH 
GO  TO  200 

TERMINATE  THE  LINE  SEARCH.  CHECK  FOR  CONVERGENCE 
400  CONTINUE 

IF(IPRN.GE.2)  CALL  PRNOUT ( 1 , NDIM ) 

IF(LEND.EQ.O)  GO  TO  440 
IF(LPAS.EQ. 1. 0R.F2.LE.F1 )  GO  TO  420 
F2=F1 
S2=S1 

DO  410  11=1 , NDIM 
R7=X1 (II) 

X1(II)=X2(II) 

X2(II)=R7 
R7=G1 (II ) 

G1 (II )sG2(II ) 

410  G2(II )«R7 
GO  TO  440 
420  F2=FM 
S2=SM 

DO  430  11*1, NDIM 
DX(II)*DM(II) 

430  X2(II)*XM(II) 

440  CONTINUE 
DF*FEND-F2 
JC0N*0 
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CALL  CONCHY (DF,S2,DX,NDIM) 

IFdPRH.EQ.  1.AND.  JCON.EQ.  1) 

•  CALL  PRN0UT(1 ,NDXM) 

IF(JCON.EQ.I)  CALL  PRN00T(5,NDZH) 

RETURN 

500  CONTINUE 
JCONsI 

CALL  PRN0UT(1,NDXM) 

IF( (NF-IX).LE.NFLX)  CALL  PRN0UT(2,NDIM) 
IF( (NF-IZ).GT.NFLX)  CALL  PRN0UT(4,NDIM) 
RETURN 
END 


SUBROUTINE  BFGSUP(S2,IH,NI,NDIM) 

THIS  SUBROUTINE  UPDATES  THE  RESSIAN  MATRIX  BY 
MEANS  OF  THE  BFGS  ( BROYDEN -FLETCHER -GOLDFARB-SH ANNO ) 
UPDATING  FORMULA.  THE  NEW  LINE-SEARCH  DIRECTION 
IS  ALSO  COMPUTED. 

WITTEN  BY  J.S.  DE  NAT ALE, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

VERSION  I:  MAY  1982. 

C0M40N/BLK5/D1(17),D2(1?),G1(17),G2(17),DG(17),DX(17)f 
•  H1(17,17),H2(17,17) 

DIMENSION  SV(17),SM(17,17) 

UPDATE  THE  INVERSE  HESSIAN  MATRIX 

DO  100  11*1 , NDIM 
100  DG(II )=G2(II )— G1 (II ) 

CALL  MVPROD(H1 ,DG,SV,NDIM) 

ALPHrDTPROD (DG.SV, NDIM ) 

AS*  1,0 

IF(ALPH.GT.O.O)  AS=-AS 
BETAxDTPROD (DX , DG , NDIM ) 

IF(IH.BQ.O)  GO  TO  120 
IF(NI.EQ.I)  IC=0 

IF(NI.GT. 1.AND.IC.EQ.0)  GO  TO  120 
GAMA*  -BETA/ALPH 
DO  110  11*1, NDIM 
DO  110  JJ*1,NDIM 
110  HI (XI, JJ)*H1 (II , JJ)*GAMA 
120  CONTINUE 

DO  130  11*1, NDIM 

1 30  SV (II ) *SQRT ( ABS ( ALPH ) )• (DX (II ) /BETA -SV (II ) /ALPH ) 

CALL  WPROD  (SV ,  SV ,  SM,  NDIM ) 

DO  140  11*1, NDIM 
DO  140  JJsl.NDIM 

140  H2(II,JJ)*H1(II,4J)-SM(II,JJ)»AS 
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CALL  MVPR0D(H1,DG,SV,NDIM) 

CALL  WPROD(SV,DG,SM,NDIM) 

CALL  M4PR0D (SM, H 1 , SV, NDIM ) 

DO  150  11=1, NDIM 
DO  150  JJsl.NDIM 

150  H2(II ,  JJ )*B2(II , JJ )-H1 (II, J J )/ALPB 
CALL  VVPFOD ( D1 , DX , SM , NDIM ) 

DO  160  II* 1, NDIM 
DO  160  JJsl.NDIM 

160  H2(II,JJ)*H2(II,JJ)-SM(II,JJ)/BBTA 

COMPUTE  A  NEW  SEARCH  DIRECTION 

CALL  MVPROD (H2 , G2 , D2 , NDIM ) 

IF  THE  NEV  COMPUTED  SEARCH  DIRECTION  IS  NOT 
A  DESCENT  DIRECTION,  SET  THE  INVERSE  HESSIAN 
TO  THE  IDENTITY  MATRIX  AND  RESTART  THE  SEARCH 

IC=0 

S2=DTPR0D (G2.D2, NDIM ) 

IF(S2.LT. 0. 0)  GO  TO  190 
IC=1 

DO  170  11=1, NDIM 
DO  170  JJsl.NDIM 
170  H2(II,JJ)=  0.0 
DO  180  11=1, NDIM 
IG=II 

G2(IG)=GFUN(IG, 1 .NDIM) 

180  H2(II,II)=-1.0 

CALL  MVPROD (H2,G2,D2, NDIM) 

CALL  PRNOUT ( 3 , NDIM ) 

S2=DTPF0D (G2,D2, NDIM ) 

190  CONTINUE 
RETURN 
END 


SUBROUTINE  CONCHK(FV,SV,XV,ND) 

THIS  SUBROUTINE  CHECKS  FOR  CONVERGENCE. 

COMMON/BLKA /EPF A , EPSA , ICON , JCON , IPRN , ICUB , XE ( 1 7 ) 
DIMENSION  XV (17) 

JC0N=1 

GO  TO  (100,200,300),  ICON 
100  DO  110  11=1, ND 

IF(ABS(XV (II ) ).GT.XE(II ) )  JC0N=0 
IF( JCON.EQ. 0)  RETURN 
110  CONTINUE 
RETURN 

200  IFUBS(FV).OT.EPFA)  JC0N=0 
RETURN 
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300  XF(ABS (SV ) .  GT . EPSA )  JCON=0 
RETURN 
END 


SUBROUTINE  PRNOUT(ZP,ND) 

THIS  SUBROUTINE  CONTROLS  THE 

PRINTING  OF  THE  MINIMIZATION  SEARCH  RESULTS. 

C0MM0N/BLK3/NI , NF, NG, PI , F2, S 1 , S2, II ( 17) , X2( 17 ) 
C0IW0N/BLE7/SFUN , SF ( 1 7 ) , XL ( 1 7 ) , XU ( 1 7 ) 

DIMENSION  XS(17) 

GO  TO  (100, 200, 300,400,500 ) ,  IP 
100  DO  110  lid  ,ND 

XS(II )=X2(II)/SF(II) 

110  XS(II)sXL(IlMXU(II)-XL(II))«(SIN(XS(II)))M2 
IF(ND.LE.6)  WRITE (6, 120)  NI,NF,NG,F1,F2,S1,S2, 

■  (XS(II),II>1,ND) 

IF(ND.GT.6)  WRITE(6, 130)  NI,NF,NG,F1 ,F2,S1 ,S2, 

•  (XS(II),IIr1,ND) 

120  F0RMAT(5X, 13,214, 1P10E11. 3) 

130  FORMAT (5X, 13,214, 1P10E11. 3/60X, 6E1 1. 3) 

RETURN 

200  VRITE(6,210) 

210  F0RMAT(/5X, '•••••  NO  ADDITIONAL  \ 

•  'FUNCTION  EVALUATIONS  ARE  PERMITTED  •••*•'/) 
RETURN 

300  WRITE(6,310) 

310  FORMAT (/5X, '•••••  THE  INVERSE  HESSIAN  HAS  ', 

•  'BEEN  RESET  TO  THE  IDENTITY  MATRIX  ••••••/) 

RETURN 

400  WRITE (6, 4 10) 

410  F0RMAT(/5X, '••••*  THE  SEARCH  HAS  BEEN  TERMINATED  ', 

•  'DUE  TO  INSUFFICIENT  PROGRESS  ••••••/) 

RETURN 

500  WRITE(6,510) 

510  F0RMAT(/5X, '•••••  THE  SEARCH  HAS  BEEN  ', 

•  'TERMINATED  AS  A  RESULT  OF  CONVERGENCE  •••••• 

RETURN 

END 


FUNCTION  DTPR0D(V1,V2,ND) 
DIMENSION  V1(17),V2(17) 
DTPRODsO.  0 
DO  100  11*1, ND 

100  DTPR0D«DTPR0D+V1 (II )»V2 (II ) 


ooooo  ooooo  o  o  o  o  o  o  o  o  o 


SUBROUTINE  VVPROD(V1,V2,V3,ND) 
DIMENSION  ?1(17),V2(17),V3C17,17) 
DO  100  11=1, ND 
DO  100  JJ=1,ND 
100  V3(II,JJ)=V1vII)»V2(JJ) 

RETURN 

END 


SUBROUTINE  MVPR0D(V3, VI ,V2,ND) 
DIMENSION  V1(17),V2(17),V3(17,17) 
DO  100  11=1, ND 
V2(II)=0.0 
DO  100  JJ=1 , ND 

100  V2(II)=V2(II)«-V3(II,JJ)#V1(JJ) 
RETURN 
END 


SUBROUTINE  MMPROD{V3,V4,V1,ND) 

DIMENSION  V1(17),V3(17,17),V4(17,17) 

DO  110  JJ=1,ND 

DO  100  II=1,ND 

V1(II)=0.0 

DO  100  KK=1,ND 

100  V1(II)=V1(II)+V3(II,IE)*VIKK,JJ) 

DO  110  11=1, ND 

110  V4(II,JJ)=T1(II) 

RETURN 

END 


FUNCTION  FFUN(ZX) 

THIS  FUNCTION  SUBROUTINE  ESTABLISHES 
THE  VALUE  OF  THE  OBJECTIVE  FUNCTION  AT 
THE  SPECIFIED  LOCATION. 

C0MM0N/BLK7/SFUN , SF ( 1 7 ) , XL ( 1 7 ) , XU ( 1 7 ) 
DIMENSION  XX(17) 

CALL  EVAL(XX,FF,0} 

FFUN*FF§SFUN 

RETURN 

END 


FUNCTION  GFUN(IV,IR,ND) 

THIS  FUNCTION  SUBROUTINE  ESTIMATES  THE  GRADIENT  OF 
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THE  OBJECTIVE  FUNCTION  BY  MEANS  OF  FINITE  DIFFERENCES. 
THE  LOGIC  FOLLOWS  THAT  OUTLINED  BY  STEWART  (1967). 

WRITTEN  BY  J.S.  DE  NAT ALE, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

VERSION  I:  MAY  1982. 

COMMON /BLK 1 /ETAF, ETAX , DMIN , FDER 

C0MM0N/BLK3/NI, NF,NG,F1 ,F2,S1 ,S2,X1 (17),X2(17) 

C«f10N/BLK5/D1(17),D2(17),G1(17),G2(17),DG(17),DX(17), 

•  H1(17,17),H2(17,17) 

DIMENSION  XX(17) 

COMPUTE  INITIAL  GRADIENT  VECTOR.  RECOMPUTE 

INITIAL  GRADIENT  VECTOR  IF  THE  SEARCH  HAS  BEEN  RESTARTED 

DO  100  11*1, ND 
100  XX(II)=X2(II) 

IF(IR.EQ.O)  GO  TO  110 

DD=DMIN 

GO  TO  160 

110  IF(NG.GT.O)  GO  TO  120 
XX(IV)*X2(IV)+DX(IV) 

ZlsFFUN(XX) 

GFUN* (Z 1-F2 ) /DX(IV ) 

NF=NF+1 

RETURN 

COMPUTE  SUBSEQUENT  GRADIENT  VECTORS 

120  CONTINUE 
R1*1.0E-10 
R2* 1 . 0E+00 

IF(ABS(G1 (IV) ).LT.R1 )  G1(IV)*R1«SIGN(R2,G1(IV)) 

ESTABLISH  GOVERNING  ERROR  BOUND 
ETAAsETAF 

Z 1 =ETAX*ABS (G 1 (IV )»X2 (IV ) /F2 ) 

IP(ETAA.LT.ZI)  ETAA*Z1 

ESTABLISH  DIFFERENCING  INTERVAL 

Z1*G1(IV)*»2 

Z2*ETAA»ABS (H2 ( IV , IV ) «F2 ) 

IF(Z1.LT.Z2)  00  TO  130 

Z1*ABS(F2/H2(IV,IV))«0.5 

D0*2.0*Z1 

Z1*D0*ABS(H2(IV,IV) ) 

Z2*  3. 0»Z  1  .  0*ABS  (G 1  (IV  )  ) 

DD*DO* ( 1 . 0-Z 1 /Z2 ) 

00  TO  1A0 

130  Z1*ABS(F2»G1(IV)/H2(IV,I?)»»2)»»(1. 0/3.0) 
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D0=2.0»Z1 

Z1=2. 0*ABS(G1 (IV) ) 

Z2= 2 . 0* ABS ( H2 ( I V , IV ) ) *D0*2 . 0*Z 1 
DD=D0* ( 1 . 0-Z 1 /Z2 ) 

140  IF(ABS(DD). LT.DMIN)  DD=DMIN»SIGN(R2,DD) 
DD=DD«SIGN ( R2 , G 1 ( IV ) ) •SIGH (R2 , -H2 ( IV , IV ) ) 

DECIDE  ON  THE  HOST 

APPROPRIATE  DIFFERENCING  FORMULA  TO  EMPLOY 

Z1*ABS(H2(IV,IV)*DD/G1 (IV))/2.0 
IF(ZI.GE.FDER)  GO  TO  150 

FORWARD  DIFFERENCES 

XX(IV)=X2(IV)+DD 

Z1=FFUN(XX) 

GFUN=(Z1-F2)/DD 

NF=NF*1 

RETURN 

CENTRAL  DIFFERENCES 

150  AA=ABS(H2(IV,IV))/2.0 
BB=ABS(G1 (IV) ) 

CC=-ETAA»ABS (F2 ) /FDER 

DD= (2. 0»CC ) / ( -BB-SQRT ( B0«BB-4 . 0«AA«CC ) ) 

IF(ABS(DD). LT.DMIN)  DD=DMIN*SIGN (R2 , DD ) 

160  XX(IV)=X2(IV)+DD 
Z1=FFUN(XX) 

NF=NF+1 

XX ( I V ) =X2 (I V ) -DD 
Z2=FFUN (XI ) 

NF=NF+1 

GFUN=(Z1-Z2)/(2.0»DD) 

RETURN 

END 


SUBROUTINE  OPEN 

THIS  SUBROUTINE  ASSIGNS  AND  OPENS  FILES  5 

AND  6.  PROMPTS  ARE  SENT  TO  THE  TERMINAL  (FILE  4). 

DIMENSION  IFIL(21 ) 

WRITE (4 ,100) 

100  F0RMAT( IX, 'INPUT  DATA  FILENAMES?' ,%) 

READ(4, 1 10)  IFIL 
110  F0RMAT(21A2) 

OPEN (UNITs5 , N  AME*IFIL , TYPE  = ' OLD ' ) 

WRITE (4, 120) 

120  FORMAT( IX, 'OUTPUT  DATA  PRINT  FILENAME*? ',$) 

READ (4, 110)  IFIL 
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OPEN (UNIT*6,NAMEsIFIL, TYPE*’ NEW1 ) 

RETURN 

END 


SUBROUTINE  EXIT 

THIS  SUBROUTINE  CLOSES  PILES  5  AND  6 
AND  TERMINATES  EXECUTION  OF  THE  PROGRAM. 

CLOSE (UNITs 5) 

CLOSE (UNITs6) 

STOP 

END 


THE  REMAINING  TWO  SUBROUTINES  —  CLAY  & 
BOUNDS  —  ARE  A  NUMERICAL  IMPLEMENTATION 
OF  THE  GOVERNING  CONSTITUTIVE  EQUATIONS 
FOR  THE  BOUNDING  SURFACE  SOIL  PLASTICITY 
MODEL 


SUBROUTINE  CL A Y ( IDIM , INC , ITNO , PROP , STOR , SIGBM , EPM , 

•  DSIGM,DEPM, C , UB.DLTAU , GAM, KIND, LARGE ) 

THIS  SUBROUTINE  EVALUATES  YANNIS  DAFALI AS ' 

BOUNDING  SURFACE  PLASTICITY  MODEL  FOR  CLAY  SOILS. 

WRITTEN  BY  L.R.  HERRMANN, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

LAST  REVISED:  JULY  1982. 

DIMENSION  PR0P(21 ),STOR(6},SIGB(6),DSIG(6),DEP(6),C(6,6), 

•  SB(3,3),SF(3,3),II(6),DLTA(3,3),DEPM(6), 

•  SIGBM (6 ) , DSIGM ( 6 ) , DEPT ( 3 , 3 ) i EPM ( 6 ) , EPB ( 6 ) 

DATA  11/11,22,33,12,13,23/.  DLTA/1.0,3#0.0, 1.0,3#0.0, 1.0/ 
ALFUN  (CV ,  RT.SINV )  »2. 0»RT*CV/  ( 1 . 0-.RT-  ( 1 . 0-RT )  «SINV ) 

SMALL  sO . 000 1 *PROP ( 8 ) 

DO  100  Id, 6 
SIGB(I)sSIGBM(I) 

DSIG(I)sDSIGMCl) 

EPB(I)cEPM(I) 

100  DEP(I)sDEPMd) 

IFdTNO.GT.  1 )  GO  TO  120 
IFdNC  .GT.1)  GO  TO  110 

INITIALIZE  HISTORY 


i 
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STOR(  1 )=PR0P(21 ) 

STOR ( 2 ) =STOR ( 1 ) 

STOR  ( 3  )=0. 5»  (SIGB  ( 1  )t-SIGB  ( 2 ) ) 
ST0R(4)=0.01*PR0P(8) 

ST0R(5)=0.0 
GO  TO  120 

UPDATE  HISTORY 

110  STOR ( 1 ) =ST0R ( 2 ) 

STOR (3 )=ST0R (3 )+ST0R (4 ) 

STOR ( 5 )*ST0R ( 5 )+ST0R (6 ) 

CONVERT  FROM  PLANE  STRAIN  TO  3-D 

120  IF(IDIM.EQ. 3)  GO  TO  140 
SIGB(4)=SIGB(3) 

SIGB(3)=STOR(3) 

DSIG(4 )=DSIG(3) 

DSIG(3)=STOR(4) 

DEP(4)=DEP(3) 

DEP(3)=0.0 

EPB(4)xEPB(3) 

EPB(3)=0. 0 
DO  130  1=5,6 
SIGB(I)=0.0 
DSIG(I)=0.0 
EPB(I)=0.0 
130  DEP(I)=0.0 

DETERMINE  3-D  INCREMENTAL  PROPERTIES 
140  GAM=PR0P(6) 

CALCULATE  EFFECTIVE  STRESS  INVARIANTS  AND 
DISTORTIONAL  STRESS  AND  CHANGE  MATRIX  COMPONENTS 
TO  TENSOR  COMPONENTS. 

XIB=0.0 

XIF=0.0 

DDILsO.O 

DILBsO.O 

DO  150  1=1,3 

DDIL=DDIL+DEP(I) 

DILB=DILB+EPB(I) 

XIB=XIBfSIGB(I) 

150  XIF=XIF*SIGB(I)«-DSIG(I) 

VOIDBs 1 . O+PROP ( 20 ) 

VOIDF«VOIDB 

IF(LARGE.EQ.O)  GO  TO  160 
VOIDB=VOIDB*EXP ( -DILB ) 

VOIDF»VOIDF*EXP ( -DILB-DDIL ) 

160  DO  170  N= 1,6 
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I*II(N)/10 

J=H0D(II(N),10) 

SB(I,J)=SIGB(N)-XIB*DLTA(I,J)/3.0 
SB(J,I)=SB(I, J) 

DEPT(I,J)*DEP(N)»(1.0*DLTA(I,J))»0.5 
DEPT(J, I)*DEPT(I, J) 

SF(I, J)=SIGB(N)+DSIG(N)-DLTA(I, J)*XIF/3. 0 
170  SF(J,I)*SF(I,J) 

GAMP=0. 0 

IF (KIND  .EQ.  0)  GO  TO  180 

GAMP*GAM 

UB=STOR ( 5 ) 

DLTAU =GAM*DDIL 
180  XIB=XIB-UB*3.  0 

XIF=XIF-(UB+DLTAU )»3. 0 

STOR(6 )=DLTAD 

SRTJBsO.O 

SRTJFsO.O 

DO  190  1=1,3 

DO  190  J=1,3 

SRT JB=SRTJB+SB ( I , J ) *SB ( I , J ) 

190  SRTJF=SRTJF+SF(I,J)«SF(I,J) 

SRTJB=SQRT ( 0. 5«SRTJB ) 

SRTJF=SQRT (0. 5»SRTJF ) 

SCOBrO.O 
SC0F=0. 0 
DO  200  1*1,3 
DO  200  J=1,3 
DO  200  K*1,3 

SCUB*SCUB+SB(I,J)»SB<J,K)»SB(K,I) 

200  SCUF=SCUF*SF(I, J)*SF(JtK)*SF(K(I) 

SC0B=SCDB/3. 0 
SC0F=SC0F/3.0 
SN3AB*0.0 

IF(SRTJB. GT. SMALL )  SN3AB* 1 . 5»SQRT(3. 0)»SCUB/SRTJB»*3 
IF(SN3AB.GT.  1.0)  SN3AB*  1.0 
IF(SN3AB.LT.-1.0)  SN3AB*-1.0 
SM3AF*0.0 

IF(SRTJF. GT. SMALL )  SN3AF* 1 . 5*SQRT(3. 0)»SC0F/SRTJF»»3 
IF(SN3AF.GT.  1.0)  SN3AF*  1.0 
IF (SN3AF. LT. - 1 . 0 )  SN3AF*-1.0 
CS3AB*SQRT ( 1 . 0-SN3AB«*2 ) 

CS3AF«SQRT( 1 . 0-5N3AF**2) 

C 

C  AVOID  ZERO  MEAN  PRESSURE 
C 

IF(ABS(XIB).OT. SMALL)  GO  TO  210 

DU*XIB 

XIB«SMALL 

IF(DU.LT.O.O)  XIB*-SMALL 
210  IF(ABS(XIF).GT. SMALL)  GO  TO  220 
DU*XIF 
XIF*SMALL 

IF(DU.LT.O.O)  XIF«-SMALL 


139 


o  o 


220  CONTINUE 


CALCULATE  ELASTIC  PROPERTIES 

DU 1 =VOIDB/3 . O/PROP ( 2 ) 

DU2s 1.5*(1. 0-2. 0BPR0P(5) )/( 1 . 0+PR0P(5) ) 

DUsXIB 

IF(DU.LT. PR0P(7) )  DUrPR0P(7) 

BB=DU1BDU 

GBsDU2*BB 

IF(PROP(5).GT.0.5)  GB=PR0P(5) 

DU 1 sVOIDF/3 • O/PROP ( 2 ) 

DU=IIF 

IF(DU.LT. PR0P(7))  DUsPR0P(7) 

BF=DU1»DU 

GF=DU2BBF 

IF(PR0P(5).GT.0.5)  GF=GB 
DO  230  M= 1 , 6 
I=II(M)/10 
J sHOD (II (M), 10) 

DO  230  N=M, 6 
K*II(N)/10 
L=HOD(II(N), 10) 

DU 1 =DLTA (K, I )*DLTA (L, J )+DLTA(K, J )*DLTA(I,L) 

C (H , N ) = (GB+GF ) BDU 1 B0. 5+0. 5B (BB+BF+2. 0«GAMP-2. 0B (GB+GF )/3. 0) 
•  *DLTA ( I , J )*DLTA (X,  L  ) 

230  C(N,M)=C(M,N) 

CALCULATE  SIZE  OF  BOUNDING  SURFACE 


XI0B=ST0R(1) 

XI0F=ST0R(2) 

XIL=PR0P(7) 

DU10= 1 . 0/(PR0P( 1 )-PR0P(2) ) 

IF (XIOB. GE. XIL . AND. XIOF. GE. XIL }  GO  TO  2*0 
XIOBSsXIOB 

IF (XIOB. LT. XIL)  XIOBSsXIL 
XIOFSsXIOF 

IF(XIOF.LT.XIL)  XIOFS«XIL 

XIOF*XIOB+DU10B0. 5#( (XIOFS*VOIDF+XIOBSBVOIDB)BDDIL- 

•  (XIOBSBVOIDB/BB+XIOFSB¥OIDF/BF)B(XIF-XIB)/3.0) 
GO  TO  250 

2* 0  XIOF=XIOB»EXP (DU10B0.5B( (VOIDBfVOIDF ) BDDIL- 

•  (V0IDB/BB+V0IDF/BF)B(XIF-XIB)/3. 0) ) 

250  ST0R(2)«XI0F 

IF ( INC+ITNO . EQ . 2 )  GO  TO  *10 


CALCULATE  BOUNDING  SURFACE  PROPERTIES 


CALL  BOUNDS ( PROP , SRT J  B , SN3AB , XSB , XIOB , XIB , GAMB , DFIB , 
•  DF JB , XKSB , DFALB , DFJ JB , BSB , VOIDB ) 

CALL  BOUNDS ( PROP , SRTJF, SN3AF, XSF,XIOF, XIF, OAMF, DFIF, 
B  DFJF , XKSF , DFALF , DFJ JF , BSF , POIDF ) 

DBsBSB-1.0 


35 


o  o  o  o 


IF(DB  .LT.  0.0)  DB=0.0 
DF=BSF-1.0 
IF(DF.LT.O.O)  DF=0.0 

CALCULATE  PLASTIC  MODULUS 

CHECK  FOR  ELASTIC  ZONE  AND  UNLOADING 

XMS=ALFUN ( PROP (17), PROP ( 1 9 ) , SN3AB ) 
DU7*0.0001«»(1.0/XMS) 

LB=0 

DDD= 1 . 0+DB» ( 1 . O-PROP (15)) 

IF(DDD.LE.O.O)  GO  TO  290 
LB=1 

HsALFUN (PROP (16), PROP ( 1 8 ) , SN3AB ) 

DUsABS (XSB) 

IF(DU.LT.DU7)  DU=DU7 
DU8=9. 0*DFIB»«2+DFJB»»2/3. 0 
DU9=II0B 

IF(XIOB.LT.XIL)  DU9=XIL 

XKB=XKSB+H»DB/DDD« ( 1 . 0>1. 0/DU**XMS )*DU8*DU9*DU10«T0IDB 

DU1=3. 0*BB*DFIB 

DU2=GB*DFJJB 

DU2P=SQRT (3-0 )*GB*DFALB 

DU3>XKB^9.0*BB*DFIB«*24GB*DFJB**24GB*(DFALB«CS3AB)**2 

SUM=0.0 

TUO.O 

IF(SRTJB.EQ.O.O)  GO  TO  280 

DO  270  1=1,3 

DO  270  J=1,3 

DUrO.O 

DO  260  1=1,3 

260  DU=DU+SB(I,K)*SB(K, J) 

T 1  *T  1 ♦  (DU- 1 . 5*SCUB*SB ( I , J ) /SRTJB**2 ) »DEPT (I, J) /SRTJB««2 
270  SUM=SUM«SB(I,J)»DEPT(I,J) 

Tl=T1-2.0»DDIL/3.0 

280  DU=  (DU  1  •DDIL-*-DU2»SUM«-DU2P*T  1  )/DU3 
IF(DU.LT.O.O)  LB=0 
290  LF=0 

DDD= 1 . 0*DF» ( 1 . O-PROP (15)) 

IF(DDD.LE.O.O)  GO  TO  330 
LF*1 

HsALFUN ( PROP ( 1 6 ) , PROP ( 1 8 ) , SN3AF ) 

DUsABS (XSF) 

IF(DU.LT.DU7)  DU=DU? 

DU8=9. 0«DFIF»»2+DFJF»»2/3. 0 
XMS= ALFUN ( PROP ( 1 7 ) , PROP ( 1 9 ) , SN3AF ) 

DU9*XI0F 

IF(XIOF.LT.XIL)  DU9=XIL 

XKF=XKSF+H*DF/DDD» ( 1 . 0*1 . 0/DU**XMS)*DU8«DU9»DU 1 0»V0IDF 

DU4=3.0«BF«DFIF 

DU5*GF»DFJJF 

DU6=XKF+9.0*BF*DFIF**2+GF,DFJF**2+GF*(DFALF*CS3AF)**2 
DU5PsGF«DFALF«SQRT (3.0) 

SUMsO.O 


136 
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11*0.0 

ZP(SRTJP.BQ.O.O)  GO  TO  320 
DO  310  1*1,3 
DO  310  J*1,3 
D0*0.0 
DO  300  K*1,3 

300  DU=DU4SF(I,K)«SF(I,J) 

T 1  *T  1  ♦  (DU- 1 . 5*SCUF*SF  (I,  J)  /SRTJP»«2  )  "DEPT  (I ,  J )  /SRTJFM2 
310  SUH=SUM+SF(I,J)«DEPT<I,J) 

T1*T1-2.0»DDIL/3.0 

320  DO* (DU4 •DDIL+DO 5  *SUM+DU  5P  #T 1 ) /D06 
IF(DO.LT.O.O)  LF=0 

CALCULATE  PLASTIC  PORTION  OF  INCREMENTAL  PROPERTIES 

330  IF(LF+LB.BQ.O)  GO  TO  410 
DO  400  M*1,6 
I»II (M)/10 
JsNOD(II(M), 10} 

DO  400  NrM,6 
K*II(N)/10 
L*MOD(II(N), 10) 

DU*0.0 

IF(LB.EQ.O)  GO  TO  360 
T2*0.0 
TlsO.O 

IF((SRTJB**4).EQ.0.0)  GO  TO  350 
DO  340  LL*1,3 
T2=T2+SB(I,LL)»SB(LL,L) 

340  T1*Tl4SB(I,LL)»SB(LL,J) 

T1=DU2P»(T1 /SRTJB«»2- 1 . 5«SCUB»SB (I, J) /SRTJB««4 -2 . 0*DLTA (I, J)/3.0) 
T2*DU2P«  (T2/SRTJB«2- 1 . 5*SCUB»SB(I,  L  )  /SRTJB«4-2. 0*DLTA  (K,L)/3*0) 
350  DU=-0.5«(DU1»DLTA(I,J)+DU2»SB(I,JHT1)»(DU1*DLTA(K,L)+ 

•  DU2»SB(E,L)+T2)/DU3 
IF(LF.EQ.O)  GO  TO  390 

360  T2*0.0 
TUO.O 

IF((SRTJF««4).EQ.0.0)  GO  TO  380 
DO  370  LL* 1 , 3 
T2*T2+SF(K, LL ) «SF (LL , L ) 

370  T1=T1fSF(I,LL)»SF(LL,J) 

T 1  *DU5P»  (  T 1  /SRTJF”2- 1 . 5»SCUF»SF  (I,  J)  /SRTJF»«4-2. 0*DLTA  ( I ,  J  )  /3 . 0  ) 
T2*DU5P»(T2/SRTJF»»2-1.5«SCUF»SF(I,L)/SRTJF»»4-2.0»DLTA(K,L)/3.0) 
380  DU*DD-0.5»(DU4»DLTA(I,J)OU5«SF(I,J)+T1)«(DU4»DLTA(K,L)+ 

•  DU5»SF(E,L)+T2)/DU6 
390  C(M,N)*DU*C(M,N) 

400  C(N,M)*C(M,N) 

410  CONTINUE 

IF(IDIM.EQ.3)  RETURN 
C 

C  CONVERT  3-D  PROPERTIES  TO  PLANE  STRAIN 
C 


DU*0.0 
DO  420  1*1,4 


DU*C(3,I)»DEP(I)4.DU 

c(3,i)*c(4,i) 

420  C(4,I)»0.0 
DO  430  1*1,3 
c(i,3)*cd,4) 

430  C(I,4)r0.0 
STOR(4)=DU 
RETURN 
END 
C 
C 
C 

SUBROUTINE  BOUNDS ( PROP , SRTJ , SN3A , X , XIO , XI , QAM , 

•  DPI, DFJ,XXS,DFAL,DFJJ,BS, VOID) 

THIS  SUBROUTINE  EVALUATES  THE  RELATIONSHIP  OF  THE 
STRESS  STATE  TO  THE  BOUNDING  SURFACE,  AS  REQUIRED  BT 
SUBROUTINE  CLAT. 

WRITTEN  BY  L.R.  HERRMANN, 

DEPARTMENT  OF  CIVIL  ENGINEERING, 

C  UNIVERSITY  OF  CALIFORNIA,  DAVIS. 

LAST  REVISED:  JULY  1982. 

DIMENSION  PR0P(21 ),FSS(3) 

ALFUN(CV,RT,SINV)*2. 0»RT»CV/( 1 . 0+RT-( 1 . 0-RT)«SINV) 
DFUN (FUN, RT, FUNC ) *FUN««2» ( 1 . 0-RT ) / ( 2. 0»RT«FUNC ) 

XN= ALFUN ( PROP ( 3 ) , PROP (4 ) , SN3A ) 

DNAL=DFUN (XN , PROP (4 ) , PROP (3 ) ) 

RrALFUN ( PROP ( 9 ) , PROP ( 1 2 ) , SN3A ) 

DRAL=DFUN (R , PROP ( 1 2 ) , PROP (9 ) ) 

A= ALFUN (PROP ( 1 0 ) , PROP ( 1 3 ) , SN3A ) 

DAALsDFUN (A, PROP( 13) , PROP( 10) ) 

YS=R*A/XN 
CC=PR0P(14) 

IF(CC.GT. 0.999)  CCsO.999 
SHIFT  PROJECTION  POINT 
D1sXI-XIO»CC 

IF(ABS(D1  ).LT. 0.001 )  DU0.001 
D2*CC-1.0/R 
D3*D1«D2 
D5*CC»(CC-2.0/R) 

Q  *SRTJ/D1 
QC* 1 . 0E+20 

IF(R*CC.NE. 1.0)  QC*XN/(1.0-R»CC) 

QO* 1 . 0E+20 

IF(CC.NE.O.O)  Q0*XN«(SQRT(1.0+IS«YS)-(1.04YS))/R/CC 
IF(SRTJ.NE.O.O)  GO  TO  100 
IF(D1.GT.  0.0)  GO  TO  120 
GO  TO  140 

100  IF(CC.LT. 1.0/R)  GO  TO  110 
IF(Q  .GE.  0.0)  GO  TO  120 


non 


IF(Q  .LB. 
IF(Q  .GE. 
GO  TO  130 
110  IF(Q  .GE. 
IF(Q  .GE. 
IF(Q  .LE. 
GO  TO  140 


QC)  GO  TO  120 
00)  GO  TO  140 

QC)  GO  TO  130 
0.0)  GO  TO  120 
00)  GO  TO  130 


PROJECTION  ON  ELLIPSE  1 

120  D4sD1«DH((R-1.0)»SRTJ/XN)«2 

BSrXIO* ( -D3+SQRT (D3#D3-D4  *(D5*(2.0-R)/R)))/D4 

LOSTd 

GO  TO  150 

PROJECTION  ON  HYPERBOLA 

130  D6=SRTJ*(1.0/R+A/XN)/XN 
D7=D3+D6 

D8=D 1 #D 1 - (SRTJ /XN )  M2 
BSa-O. 5*XIO* (D5-2. 0«A/R/XN )/D7 
IF(D8.EQ.0.0)  GO  TO  150 

BS=XIO»(-D7+SQRT(D7»D7-D8«(D5-2. 0«A/R/XN) ) )/D8 

L0STr2 

GO  TO  150 

PROJECTION  ON  ELLIPSE  2 

140  T=PROP( 1 1 ) 

FOP=XN/SQRT( 1 . 0+TS»*2 ) 

XJO=A* ( 1 . 0+YS-SQRT ( 1 . 0+YS»»2 ) ) /YS 
BT*T« (XJO-T»FOP ) / (IJO-2. 0»T»FOP ) 

RO= (BT-T ) /FOP/XJO 
PSIal . 0/(R*(BT-T) ) 

D9=T-BT+CC 

D10=D1*D9 

D 1 1 aD 1 «D 1 ♦RO*SRTJ •SRTJ 

BSaXIO* (  -D 1 0+SQRT (D 1 0*D 1 0-D 1 1»(D9*D9-BT»BT) ) )/D1 1 
LOSTa3 

150  XIBARaBS»(XI-XIO»CC)*XIO»CC 
IF(XIBAR.EQ.O.O)  XIBARal. OB-20 
TH=BS*SRTJ /XIBAR 
XaTH/XN 
DOaXIO 

IF(XIO.LT.PROP(7))  DU»PROP(7) 

D0Sa12. 0»VOID/(PROP{ 1 )-PR0P(2) )«XI0«»2»D0 
GO  TO  (200,300,400),  LOST 

NORMAL  CONSOLIDATION  ZONE 


200  PSIaYS/(R-1.0)«*2 

DU*R*  ( 1 . 0+XUfR*  (R-2.0)*X*X) 

GAMa ( 1 . 0+ (R- 1 . 0 )»SQRT ( 1 . 0*R» (R-2 . 0 )«X«X ) ) /DO 
DFI*2. 0»XI0» (GAM- 1 . 0/R  )»PS1 


o  n  n  n  o  o  o  o  o  o 


3 


DFJJ=2. 0*XIO*GAM* ( (R-1 . 0)/XII  )**2*PSI*BS/XIBAR 
DFJ  *DFJ J •SHTJ 

XKS=D0S*(GAM-1. 0/R)*(GAM+R-2. 0)*PSI*PSI/R 
DFAL*PSI*6 . 0* (R- 1 . 0 ) *TH*GAM*XIO* ( ( (R- 1 . 0)/ (R**2* 

•  (2. 0/R-GAM-1 . 0)  M .  0  )*DRAL-(R-1 . 0)*DNAL/XN)/XN**2 
RETORN 

OTERCONSOLIDATION  ZONE 
300  DUs 1 . 0-X* ( 1 . O+YS ) 

GAH=- (DU+SQRT ( ( X-YS- 1 . 0 )**2* ( X*X- 1.0) *YS*YS ) ) / (R* (X*X- 1.0)) 
DFI=2. 0*XIO* (GAM-1 . 0/R ) 

DFJ=2. 0*XI0*(( 1 . O+YS ) /R-X*GAM ) /XN 
XKS=DUS* (GAM-1. 0/R)*(DU*GAM+2.0*A/XN)/R 
DFJJsDFJ/SRTJ 

DFAL=6. 0*XIO*(DNAL*(TH*GAM/XN-1. 0/R+A/(R*TH*GAM)-2. 0*A/XN )/ 

•  XH**2+DRAL*( 1 . O/TH-1 . O/XN+A/ (XN*TB*GAM ) ) /R**2+DAAL* ( 

•  1 . 0/XN-1 . 0/(TH*GAM*R ) )/XN ) 

RETURN 


TENSION  ZONE 

400  GAM= ( -T+BT-SQRT (BT*BT-RO*TH*TH*T* (T-2 . 0*BT ) ) ) / ( 1 . 0+R0*TH*TH ) 
DFI=2. 0*PSI*XI0* (GAM+T-BT ) 

DFJJ=2. 0*PSI*XXO*GAM*RO*BS/XIBAR 
DFJ=DFJJ*SRTJ 

XKS=DUS*PSI*PSI* (GAM+T-BT )• (GAM*  (BT-T )+T* (2. 0*BT-T ) ) 
DYSAL=YS* (DRAL/R+DAAL/A-DNAL/XN ) 
DFOPAL=FOP*(DNAL/XN-YS*DYSAL/( 1 . 0+IS*YS) ) 

DJOALrXJO* (DAAL/A-DYSAL/YS )+A* ( 1 . O/YS-FOP/XN ) *DYSAL 
DBTAL= ( (T-BT ) *DJOAL- (T-2. 0*BT )*T*DFOPAL ) / (XJO-2 . 0*T*P0P ) 
DROAL=DBTAL/FOP/XJO-RO* (DFOPAL/FOP+DJOAL/XJO ) 

DFAL=3. 0*PSI*XI0*TH*GAM*  (DROAL+2 . 0*RO*DBTAL/ (T+GAM*  2 . 0*BT ) ) 

RETURN 

END 


c 

c 


END  OF  MODCAL 


i 
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